PERTURBATION  ANALYSIS  OF  NONLINEAR 
REACTOR  KINETICS  WITH  DOPPLER  FEEDBACK 


Edward  Paul  Nemeth 


B%8  \ 


hool 


Monterey,  California 


T 


**  I 


^  ■■■-'    Li 


Perturbation  Analysis  of  Nonlinear 
Reactor  Kinetics  with  Doppler  Feedback 


by 


Edward  Paul  Nemeth 


Advisor: 


D.  PI .  Nguyen 


June  1972 


Approved  {^oh.  public  k.(lIqx&<l',    dU&ubvutlon  unLimitcd. 


Perturbation  Analysis  of  Nonlinear 
Reactor  Kinetics  with  Doppler  Feedback. 


by 


Edward  Paul  Nemeth 

Ensign,  United  States  Navy 

B.S.,  Villanova  University,   1971 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  MECHANICAL  ENGINEERING 


from  the 

NAVAL  POSTGRADUATE  SCHOOL 
June  1972 


c-  ^ 


ABSTRACT 

The  analysis  of  nonlinear  reactor  kinetics  with  Doppler  feedback 

is  developed  by  the  use  of  the  Poincare-Lindstedt  perturbation  technique. 

The  two  types  of  temperature  dependent  reactivity  coefficients  to  be 

investigated  are 

dk      ^        I  ,,  dk  1 

dT      *         T  6nd  dT     "         T   3/2    . 

Each  of  these  reactivity  models  are  analyzed  for  the  reactor  with  positive 
and  negative  coefficients  of  reactivity  and  infinite  slab  and  finite  cylin- 
drical geometries.    The  Poincare-Lindstedt  perturbation  solution  is  con- 
structed in  the  form  of  an  expansion  about  the  unperturbed  state  of  the 
reactor.    Each  of  the  perturbed  equilibrium  solutions  is     subsequently 
analyzed  for  stability. 
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I.     THE  NONLINEAR  TEMPERATURE  DEPENDENT  REACTOR 

A.         GOVERNING  EQUATIONS 

In  a  homogeneous,  one-velocity,  bare  reactor  the  neutron  popu- 
lation is  described  by  the  equation  (1) 

±  m^i)  =  p  vv*,-o  -  ia  m,-i)       (i) 

i 

The  first  term  Eq.   (1)  represents  the  time  rate  of  change  of  the  neutron 

flux  within  a  differential  control  of  volume.    This  term  is  balanced  by  the 
net  effect  of  the  diffusion  of  neutrons  ,  the  absorption  of  neutrons  ,  and 
the  production  of  both  prompt  and  delayed  neutrons  within  the  volume. 
The  concentration  of  the  precursors  is  described  by  [1] 

i    =    1,2,3,  ■■■    tN        . 

The  solution  of  these  equations  with  temperature  dependent  parameters 
represents  the  goal  of  this  work. 

The  temporal  and  spacial  coordinates  may  be  nondimensionalized 
by  introducing  the  definitions 


-  V  Vd 


x   =  v  *y  d   #    and       r  =  at  f*  -t   . 


Equations  (1)  and  (2)  become  the  following 

||  =  v>  -  ip  +  i  (L-p)  <f  +  %-£_  x  Cl  <3) 

where  the  precursor  concentration  and  the  flux  must  have  a  zero  value 
on  the  boundaries  of  the  reactor.    The  solution  is  desired  for  the  infinite 
slab  reactor  and  the  finite  cylindrical  reactor.    The  first  of  these  cases 
is  the  slab  which  may  be  envisioned  as  a  reactor  which  is  infinite  in 
two  dimensions  and  bounded  by  the  coordinates  +1  and  -1  in  the  third 
dimension.    This  simple  geometry  will  easily  lend  itself  to  analysis, 
but  permits  the  observance  of  mathematical  rigor  and  the  physical  inter- 
pretation of  the  results.     The  boundary  conditions  require  that  the  flux 
is  zero  at  the  coordinates  +  1.    The  second  case  is  that  described  by 
a  finite  cylinder  which  has  boundaries  on  the  longitudinal  axis  at  +  1 

and  at  a  nondimensional  radial  distance  of  r    .    Therefore,  the  flux  is 

o 

required  to  be  zero  at  the  extremes  of  the  longitudinal  coordinate, 

Z  =  +  1 ,  and  at  the  radius  r  =  r    . 

~  o 

B.         DOPPLER  EFFECT 

The  lumped  parameter    k    is  generally  referred  to  as  the  neutron 
multiplication  factor.    In  elementary  reactor  analysis  the  multiplication 
factor  is  evaluated  as  a  constant  given  by  the  product    [l] 


*      -      *.   Pr* 


However,  in  reactors  that  operate  at  high  temperatures  the  multiplication 
factor  becomes  temperature  dependent  as  will  subsequently  be  shown. 
One  of  the  major  considerations  to  be  made  when  accounting  for 
the  temperature  dependence  of  the  multiplication  factor  is  that  described 
by  the  nuclear  Doppler  effect.     Physically,  this  phenomenon  results  from 
the  temperature  broadening  of  the  capture  cross-sections  of  the  fuel 
nucleii.    Thus  the  Doppler  effect  may  be  envisioned  by  a  temperature 
dependent  competition  for  neutrons  by  capture  and  fission  processes. 
The  quantitative  description  of  the  Doppler  effect  is  achieved  by  making 
the  following  assumptions: 

a.  The  fuel  density  as  a  function  of  velocity  is  given  by  a 
Maxwellian  distribution  function. 

b.  The  capture  cross-sections  near  an  isolated  resonance  is 
given  by  the  Breit-Wigner  formula . 

It  may  be  shown  that  the  average  radiative  cross -section,      O"  -g        , 
and  scattering  cross-section,     0*s         /  are  described  as  a  function  of 
energy  and  temperature  given  by  [1] 

cra  (£J)      =     -!—*    f(i  A) 

and 


where       j    (    5  >  #)     and      a-  C  5>X/       are  tabulated  functions  of 
temperature,  o 

r 


v       A 

The  effect  of       0"y       and      (T  $         is  most  conveniently  incorporated 

into  one  quantity  by  introducing  the  effective  resonance  integral,  I, 
defined  by  [2] 


j     =     fC(T)  dg        f     oyo; 

4-  j  ^a^rcECT.vt    e       J  a;  +  cr5  + 


—  ^1 


where       O"^       is  the  additional  scattering  cross-section  per  fuel  atom. 
It  is  usually  more  convenient  to  analyze  the  temperature  dependence 
through  the  J  function  which  is  related  to  the  resonance  integral  as 
follows  [1] 


I    =     %r*     J(hfi) 


where  M*,  ov„   V 


f3  NA  cr,    ty 


and  the  subscript  M  refers  to  the  moderator. 

If  attention  is  returned  to  the  concept  of  the  multiplication  factor, 
its  dependence  upon  temperature  may  be  investigated  [1] 
Since 


•f 


10 


and 

jo    «      e^p  (  I)      > 

the  Doppler  effect  on  k  is  studied  by  considering  the  temperature  effect 
on  JO     .    However,  the  preceding  analysis  showed  that 

i  -    n?*  SW 

and  it  may  be  simply  concluded  that 


\  «    e*p   fja>/3)    • 


If  the  derivative  with  respect  to  temperature  is  taken,  it  is  obvious  that 

However,  under  most  operating  conditions  the  value  of     ^     is  near 
unity.    Thus  for  small  changes  in      R.     the  following  relation  is  valid, 

di       ^         dJ  t-\ 
oc •  (5a) 

dT  dT 

The  left  hand  side  of  Eq.    (5)  is  referred  to  as  the  temperature  coefficient 
of  reactivity  and  provides  a  functional  dependence  upon  the     J       function 
and,  in  turn,  upon  the  changes  in  temperature  of  the  reactor. 

In  a  paper  by  A.  Reichel  the  temperature  dependence  of  the 
function  was  investigated  [2].     The  results  were  based  upon  evidence 

o  o  o 

gathered  from  experimentation  with  U  and  from  analytic  investigations 


11 


into  the  nature  of  the     <J      function.    Reichel  concluded  that  physically 

realizable  states  of  the  reactor  must  be  located  in  the  range 
.02   c     $     C   .6       and  .  OOOi    <  ft  <    10   • 

The  results  are  most  simply  presented  by  Figure  1.    The  variable 
/3        is  described  by     K       where  K     =  V /- 

The  shaded  region  represents  states  which  are  physically  realizable  for 

900 

U        .    From  the  analytic  investigations  it  was  observed  that  values  of 

1$      at  the  extremes  of  its  range  coincided  with  conditions  which 
required  that  the      0       function  be  a  constant.    Therefore,  the  tempera- 
ture dependence  of      J        is  destroyed  and,  in  turn,  the  Doppler  effect 
and  temperature  dependence  of  the  coefficient  of  reactivity  is  nonexistant. 
In  other  regions  of  Figure  1  the  dependence  of      J       was  found  to  vary 

with  temperature  raised  to  a  constant  power.    These  regions  are  indicated 

const 
on  Figure  1  by  the  lines  labeled  with  T  .    Thus  for  any  particular  state 

of  operation  defined  by  values  of     AS      and       5       the  dependence  of 
upon  temperature  may  be  defined.     If  Eq.   (6)  is  evoked,  this     J       de- 
pendence is  transferred  to  the  governing  equations,  Eqs.   (3)  and   (4), 
via  the  temperature  coefficient  of  reactivity. 

As  may  be  noted  from  Figure  1  the  variance  of  the      J        function 
with  temperature  may  take  on  many  different  forms.    However,  in  an 
analysis  presented  by  Hummel  and  Okrent  [3]  the      0       function  was 
found  to  vary  between  the  functions 

J     «■    ^  T  J   <x     T _l/i 

12 


J(hfi) 


FUNCTION  TABLE 


-3/2 

T  Model 


Figure  1 
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These  results  are  used  to  determine  the  functional  dependence  of 
with  temperature.    The  results  are  that 


a 


and 


dT 

dT 


oc 


oc 


T 


-t 


T 


-3 


(6) 


h 


1 


respectively.       In  the  remainder  of  this  paper  the  effect  of  these  types 
of  temperature  feedback  are  analyzed. 

The  proportionalities  in  Eq.    (6)  are  removed  by  the  addition  of  a 

are  determined  by  experiments  per- 


constant.    The  values  of 


dT 


formed  on  the  reactor  in  question.     Some  representative  values  are  given 
for  various  reactors  . 


TABLE  1    [4] 


Reactor 


Nuetron 
Lifetime 


Temperature 


Coefficient 


Fermi 

EBR  II 

Hallam 

Peach  Bottom 

EGCR 

EBOR 


14  JU  sec 

68  u.  sec 

37  m   sec 

23  wi  sec 

4  0  *v\  sec 

22  jul   sec 


-3.4  x  10 
-4.2  x  10" 
+  .3x10 
-2.6  x  10 
-4.5  x  10 
-1.5  x  10* 


-5 


-5 
-5 


-5 


A  k/k  per  °C 
A  k/k  per  °C 
A  k/k  per  °F 
A  k/k  per  °C 
A  k/k  per  °C 
A  k/k  per  °C 


In  the  subsequent  analysis  the  need  to  differentiate  between  the 
feedback  models  will  arise.    The  distinction  will  be  made  by  referring 
to  the  "T_1  model"  or  "T"3/2  model." 
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-1  -3/2 

For  the  T       and  the  T  model  reactors  the  temperature  dependence  of 

the  multiplication  factor  may  be  described  by  the  equalities 

dji  a' 

dT     =     "    r 

and      41  £ 

CH  T        ~  J  Vk 

respectively.    It  may  be  noted  that  the  negative  sign  is  incorporated 

into  above  equation.  This  convention  is  generally  observed  because  of 
the  greater  frequency  of  reactors  with  negative  temperature  coefficients. 

d&/. 

In  Eq.    (7)  the  units  of      G.     and      Q.      are  defined  such  that  /  QT 

has  the  units  of       &    k/k  per  °C. 


C.        ASSUMPTIONS 

Before  the  analysis  of  the  reactor  may  be  undertaken  two  further 
assumptions  must  be  made  concerning  the  dependence  of  the  neutron 
flux  with  temperature  and  the  delayed  neutron  concentration  with  time. 
The  first  assumption  states  that  the  reactor  temperature  rises  instantan- 
eously with  the  flux  (this  is  the  usual  prompt  feedback  model)  .    Equation 
(7)  becomes 

df 5l 

and  (8) 

The  second  assumption  states  that  the  delayed  neutron  concentration 
has  delayed  neutron  concentration  has  reached  its  steady-state  value. 


15 


TABLE  2    [1] 

Delayed  Neutron 

Group 

Half  life  (sec) 

n  87 
Br 

54.5 

137 

24.4 

u   88 

Br 

15.3 

138 

6.3 

n  89 
Br 

4.4 

ou93,94 

The  justification  of  this  assumption  lies  in  the  comparison  of  the  neutron 
life  time  in  Table  1  with  the  precursor  half-life  in  Table  2  .     It  is  evident 
that  a  number  of  neutron  generations  will  have  progressed  before  the  pre- 
cursor concentration  has  been  changed. 

The  assumption  of  steady-state  precursor  concentration  reduces 
the  neutronic  eguation  to  the  case  without  delayed  neutrons.     This  work, 
therefore,  deals  with  the  case  where  either  no  delayed  neutrons  exist 
in  the  reactor  or  the  delayed  neutron  concentration  instantaneously  follows 
the  flux  distribution.    The  latter  case  is  realistic  when  considering  the 
reactor  response  to  small  perturbations. 


16 


II.     ANALYSIS  OF  THE  REACTOR  WITH  A 
TEMPERATURE  COEFFICIENT 


!~     =     -  J        (T    l  Model) 
dT  T 


The  solution  of  the  nonlinear  equations,  Eqs .    (3),   (4)  and  (8), 
governing  the  reactor  kinetics  is  accomplished  by  means  of  the  Poincare- 
Lindstedt  perturbation  technique  [5].    This  method  offers  a  solution  to 
the  neutron  flux  distribution  for  small  perturbations  in  the  flux  intro- 
duced by  the  temperature  feedback  effect.    The  perturbed  flux  distribution, 
referred  to  as        t(Kj  O        is  described  by  a  Taylor  series  expansion 
about  the  initial  state       To  (  "X)     where      k.       is  the  amount  of  perturbation 
from      7©  (*)       •     The  perturbation  technique  yields  multiple  equilib- 
rium solutions  for         tY^O         and  associated  eigenvalues  for  each 
of  these  states.    The  presence  of  more  than  one  equilibrium  state,  known 
as  the  bifurcation  phenomenon,  requires  the  investigation  of  the  stability 
of  each  of  the  possible  states  of  flux.     Stability  considerations  are  made 
by  imposing  an  initial  deviation  on  an  equilibrium  state  and  examining 
the  resulting  temporal  behavior  of  this  state. 

A.         EQUILIBRIUM  STATES  FOR  A  SLAB  REACTOR 

Before  the  general  solution  to  Eqs.    (3),    (4)  and  (8)  is  undertaken 
the  equilibrium  states  for  the  slab  reactor  is  investigated.    For  the  case 
where  the  flux  does  not  vary  with  time  (i.e.     -v  sr     —    <J     )  and  the 

17 


reactor  kinetics  are  not  dependent  upon  delayed  neutrons  the  governing 
equations  reduce  to  the  following  boundary-value  problem  in  the  spacial 
domain  D[-l ,  1]: 


<Pcx)     =    O 


>n     <SO 


and 


d  fe(tp) 
d  9 


PL 


(8) 


Equation  (8)  may  be  solved  to  give 


&(*)    -     1    -    a  JU. 


¥o<X> 


(9) 


This  result  may  be  substituted  into  the  boundary-value  problem  and  the 
result  is 


V'tftx)  +    f>  -  a  J^ 


C£o0 


where 


<P(X)     ■=     O 


A      -        ^.    -    1 


<^(x)   =    O       in   D 


(10) 
on    &  D 


The  solution  to  Eq .   (10)  is  constructed  about  a  flux  distribution 

To  (■  *  )  and  eigenvalue       A  o  which  will  be  determined. 

Physically,         A  0        is  this  eigenvalue  corresponding  to  the  initial 

material  buckling  of  the  reactor  before  perturbation  by  feedback.    After 

18 


the  system  is  perturbed  the  flux  and  the  material  buckling  establish  a 
new  state  which  may  be  described  by  a  Taylor  series  expansion  about  the 
initial  state.     If  the  amount  of  perturbation  is  measured  by      €       ,  the 
resultant  reactor  state  is 

"      „o  (U) 


Thus  for  small  perturbations  of  the  system  the  values  of  i    (  0< ,  O 

and      A(0       may  adequately  be  described  by  truncating  the  series  aft 
the  third  term.     The  unperturbed  flux,         So  ("X)        ,  and  eigenvalues, 
A  ©  /  are  defined  as  the  solution  to  the  linear  problem 


er 


V1(foM  +  (X0-  a)  <£cx)    =  O 


m   P 

(12) 
<fo(*X)       =      C  on     <5P. 


Since         A©        and      Q.      are  constants, the  solution  is 


%<*)    =  /L     Cos  V^°~ a'   * 

and  2  (13) 

A.  -  a    -     r%  ) 

where  yy\  =  1 ,  3  ,  5  ,    .  .  .    .    For  the  purpose  of  simplifying  the  subsequent 
analysis  the  shape  function        T /n       ^or  t^ie  s^a^  reactor  is  defined  as 

%        -     COS      f"&      *)  •  U4. 
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The  remaining  terms  of  the  expansion  given  by  Eq.   (11)  are  determined 
by  differentiating  Eq.   (10)  with  respect  to     £        and  taking  the  limit  as 

G        tends  to  zero.    The  values  of      y*      and      A      are  then  determined. 
By  repeating  this  process  the  values  of      *r      and      A       are  evaluated. 
Thus,  differentiating  Eq.    (10)  once  yields 

V2<£  +  if  +  \$-  a  <f  l*[  %,  -  a  ¥  =  o.  as) 

The  limit  as      C       tends  to  zero  requires  that 
As      £        tends  to  zero,  Eq.    (15)  reduces  to 

S/1(f(*)   ■+  (A„-  a)  </Vx)  =  -Af.oo      i«  a>  (i6) 

<¥>{X)~0  OKI         SD    . 

If  Eq.    (16)  is  to  have  a  nontrivial  solution,  the  solvability  condition 
(provided  in  Appendix  A)  may  be  evoked.     Thus 


/( 


-A  £00)  P"cx)    d-x     -  o 


D 

-     V 
where  ^P    C-X)  is  the  homogeneous  solution  to  Eq.    (16).     The 

solution  for         y^  (x)  is  easily  seen  to  be        7P  (x)    ~    &/n      // 


/-I 
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When  the  functions  Sr  CX)        and  To  ^)         are  substituted 

into  the  solvability  condition  the  result  is 

-A   f  /U  ^  <*    **    -  °    • 

Since  the  integral  is  nonzero,  it  follows  that  A  =  O     .    This  result 

is  substituted  into  Eq.    (16)  and  the  solution  for  ^P  (  *)  is  clearly 

seen  to  be        r(^)  =  O/^     'Ai     .    Thus  the  linear  terms  of  the  expansions 

are 

lf(x)      =        B«      ^    (X)      and 

*    =    o     . 


The  parameter      £       is  chosen  to  represent  the  average  magnitude  of 
the  linear  term  in  the  expansion  of  Eq.   (11),  weighted  by  the  i^x 

distribution,  and  measured  in  units  of  the  maximum  initial  flux.     Thus, 


Differentiating  Eq.    (17)  with  respect  to       £      yields 

A*    =  /  %W     fit*,*)     dx 


(17) 


where 


9^*,0    =      ^f*)    +    €    <£(x)    + 
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By  letting      c       approach  zero 


thus     /"1/vi     =      ^/>t     •     The  second  terms  of  the  expansions  are 

<f  =      /l^      ^  and  .     A      -    O     .  (18) 

The  third  term  of  ^P(Xi^)       and        A  ^ )        are  similarly 

determined  by  differentiating  Eq.   (16)  with  respect  to      €       .     The  result 


is 


V^Ooo    +  A  CO  </>(x,o     +  2M^>  <P (*,<:) 


ifi  C  K  £ )      -      O  on     5  »  . 

Since         A~  O    /  the  result  of  taking  the  limit  as      £       tends  to  zero 


is 


V  ¥(*)   +   (Ac-  a)  v  w    =  (is) 


?. 


/ 


<~P  CX)         ~      O  on       !>*>. 


The  solvability  condition  requires  that 


/  [a<foo  ~^T  -  \%m)  <Pi*i  dx  =  o 


3) 
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where        *jr    (x)  is  the  homogeneous  solution  to  Eq.   (19).    This 

solution  is  obviously  seen  to  be         y^  (  "X )    ~~       L  *  a\ 

When  the  various  functions  are  substituted  into  the  solvability  condition 
the  following  becomes  evident 


X 


and       /\    ~    CL         .    This  result  is  replaced  in  Eq.    (19)  and  it  is  found 
that  the  right  hand  side  identically  becomes  zero.    Therefore,  the  solution 
for  the  third  terms  of  the  expansions  is        j^CxP   —    C/y\     i/h     and        (20) 

]\    ~     CX        •    The  equilibrium  solution  for       ^P("X,£)         is  obtained 
by  placing  the  results  of  Eqs  .   (13),   (18),  and   (20)  into  the  series  of 
Eq.    (11).     Thus, 

?Yx,o   =  /)„  (l  +  £  +  i  €2  £  +<••)  %  M  (2D 

and 

Xco-  a  =  C^/2)2   +  2  £2a+-- 

where /V\  =1,  3,   5   ...    .     The  value  of        *r(X,  O        fromEq.   (21)  may 
be  evaluated  if  the  constants  A     and  C,  are  defined.     Since  the  reactor 
is  assumed  to  operate  at  its  fundamental  mode  before  temperature  feed- 
back becomes  effective  (later  it  will  be  shown  that  the  other  modes  are 
unstable) ,  the  value  of  A     may  be  determined  by  the  power  requirements 
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of  the  reactor.    Thus, 


p--u 


D 

where    P    is  the  linear  reactor  power  level  and      Q     is  a  conversion 

factor.    The  value     L/^       is  undetermined  and  will  remain  an  unknown 

l 

of  order      €  C^  throughout  the  analysis  .    For  each  mode,   /Yl       ,    Y<x,6.) 

as  defined  in  Eq.   (21),  describes  two  possible  states  in  the  reactor 
as  depicted  in  Figure  2  for  the  fundamental  mode.    The  state  of  larger 
flux  is  determined  by  a  positive  value  of      c       while  the  smaller  state 
is  described  by  a  negative  value  of      £        .    Associated  with  each  state 
of  the  reactor  is  an  eigenvalue,     A^'  Q-     ,  given  by  Eq.    (21)  and 

plotted  against       c       in  Figure  3.    The  values  of      A^£  )        CX.      for  i\\e 
case  of  negative  reactivity  are  drawn  as  a  solid  parabola  which  opens 
to  the  right.     It  may  then  be  noted  that  the  equilibrium  state  may  only 
exist  for  values  of     A(w  "     &■      which  are  greater  than  the  value  of 

\Q        ,  the  initial  material  buckling  of  the  unperturbed  system.    For 
the  reactor  with  positive  reactivity  the  converse  is  apparent.    The  cor- 
responding curve  for  positive  reactivity  in  Figure  3  is  that  drawn  with 
a  dashed  line.    These  eigenvalues  are  found  to  be  always  less  than  the 
value  of        Ao        .It  may  also  be  noted  that  the  eigenvalues  of  the 
system  are     X(0  ~   Q-      where  the  term  (     —  Q.      )  describes  the 
amount  the  material  buckling  changes  so  as  to  accommodate  the  temper- 
ature feedback. 
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FLUX  DISTRIBUTION  FOR  SLAB  REACTOR 


^P(*>£)    positive       £ 
^r(*X,0    negative     £ 


Figure  2 
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negative  reactivity 
positive  reactivity 
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B.         STABILITY  ANALYSIS  OF  THE  SLAB  REACTOR 

The  stability  criterion  for  a  slab  reactor  is  developed  by  considering 
the  temporal  effect  of  an  initial  deviation  on        r  (a;  £  )       •    Under 
these  conditions  the  governing  equations,  Eqs.    (3),    (4)  and  (9)  reduce 


to 


UW     I) 


in  1> 


(9) 


V(^O)      =      ^oCx) 


;     ^Cx,?)  =  o 


or\ 


&* 


These  equations  may  be  compacted  into  the  more  convenient  form 

in        £       (22) 
(-PiX)*)      -    O  on       iD 


where 


A    -    8.  -  1   . 


The  time  dependent  nature  of  the  neutron  flux  is  investigated 
by  prescribing  an  arbitrary  initial  deviation,    "YI        ,  on  the  equilibrium 
flux  and  then  observing  the  growth  of  the. deviation.    An  unstable  state 
of  flux  is  that  which  will  result  in  an  unbounded  growth  of  the  deviation 
A  stable  system  will  behave  in  a  neutrally  stable  manner  which  implies 
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that  the  deviation  tends  to  a  new  equilibrium  state  or  in  an  asymptotically 
stable  manner  which  implies  that  the  deviation  decays  to  zero.    The  time 
dependent  flux  is  described  by  a  Taylor  series  expansion  about  the 
equilibrium  state        ^(X;0       ,  defined  by  Eq.    (21), 

,      %  (23) 

such  that 

Equation  (22)  is  differentiated  with  respect  to      £      .    Therefore, 


it 


S   =v\+(^)-a-a.^^>,     -  » 


^     =     O  ovx       5  J) 


where 


The  limit  as     "^7      tends  to  zero  is  taken  and  the  result  is 

— ^E  =  V  %^t%)    +(^a)  -  a    (24) 
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If  a  solution  for       r^  ( X,  6,  T)       is  chosen  to  be 

9?  0,*,t;     =     u.<:x,o     e^pC^cotO  (25) 

then  the  stability  of  the  system  is  analyzed  by  inquiring  whether  the 
value  of      OL  (£)        is  negative  or  positive.    A  form  for  oL  C  O 

is  given  by  the  series 


cx^CO     =      dQ   ■+     £  ^    +    2  £*  ^   +  "  '   .    (26) 

The  value  of      OC  (  w       is  determined  by  the  perturbation  tech- 
nique previously  used.     In  this  vein,  Eq .    (2  5)  is  substituted  into  Eq.    (24) 
and  the  resulting  equation  is 

y^uc-x, c-)    +(>(*)-  a  -   c^co    -         (27) 

a       J^l  ~%7x)    J      UCX;t)     =     O  m      J> 

u  (x, e)     =     O  on    SD. 

Setting       c       equal  to  zero 

1  /  x  \  (28) 

7    U  (x,e)      -  {  Ao/y,  -    Ol   -  ^oj    UCxyo)  ■=   O       m    D 


U(XyO)    «      O 


Or\ 


<Sx> 


Since  the  values  of    Ao/m     an<^  (  Q»  )  have  been  prescribed  by  Eq.    (13), 
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Eq.    (2  8)  may  have  a  nontrivial  solution  only  if  (proof  is  given  in 
Appendix  C) 

Thus  the  value  of       0io        is  defined  as 

oLor*f»\      ~         (  Ao/*  "    CLj    -    (   ^o/wv  ""    Ow)  (29) 

=  (y2r  u>-^) 


for  odd  intergers  of    /Y[       and     /W^    where    /Wl      is  the  mode  of  the 
deviation.     If  the  fundamental  mode  of  flux  distribution,    /Vt  —  1    ,  is 
considered,  it  is  evident  that  Eq.    (2  9)  reduces  to 

d0,^    =   (7/2)L   (1  -  wl) 

Thus  the  values  of  0(.o\  /m         are 

-     ~  2TJ  />n    *  '3,5,7,  -••  . 
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For  values  of      c        near  zero  the  first  term  of  the  series  for         C*CO 
indicates  that  the  fundamental  mode  is  not  unstable  for  any  value    /Yl\ 
However,  the  neutral  stability  predicted  by    /V>1  -  J       obliges  an  investi- 
gation into  the  higher  order  terms  of  Eq.    (26).     Some  conclusive  results 
may  be  obtained  by  considering  the  linear  term  of     qL  (  £  ) 

For  higher  modes  of  the  reactor  (i.e. ,  '/Y\  —  1     )  Eq.    (2  9)  indicates 
that 

O  /m     —      Aa 


2^ 


/WV         <         /H 


Since  a  positive  value  of     (A  (O       may  occur,  it  is  concluded  that  the 
higher  modes  of  the  reactor  are  unstable.    For  a  physically  realistic 
situation,  however,  only  the  fundamental  mode  exists  in  the  reactor 
before  the  perturbation. 

The  solution  for      LX  (  'X,  O)       is  obtained  by  substituting  the 
value  of  Eq.    (2  9)  into  Eq .    (2  8).     The  result  is  that 

VZUCX,0)    +    fAo/v^-Cx)    LUX,0)    =  o        \n    0(3  0) 

UU,0)       -     O  °»      *»• 

It  is  then  obvious  that 

uc-x.o)    -    U^    '^ 
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where  the  index    /7YJ      designates  the  modes  of  the  deviation.     The  value 
/Wl  a  1       refers  to  the  fundamental  mode  of  the  deviation  while     /yy\  >  1 
refers  to  the  higher  modes. 

The  linear  term  of  the  expansion  for      0(  CO       ,  defined  in 
Eq.    (2  6),  is  derived  by  differentiating  Eq.    (2  7)  with  respect  to      t 
thus 

V1  ucx,o   +  (A CO  -  a  -  oico 

-a  hv\     "Wx)J   U(x,0    +    (  MO  -  oiU) 


on 


© 


Ul(X,0    =    o  on      S£>  . 

The  limit  as      £        tends  to  zero  is  taken  and 

«l   '  /  -x  y  \     ■  (32) 

V      UCX,0)       +     ^    Ao/m    ~    a    -      0Co^jm)     U.  CX.O)   = 

p  A(°)     +      ocCo)   +    a       <^oc;J  U(x,o)      'm   J> 

U.  CX;0)         =       O  OA      i». 

i 

By  noting  that       ]\     "=    O         the  solvability  condition  for  Eq  .   (32) 
requires  that 

f  /  •/  ~        ^^   \  '  H  * 

J  (oL  ■+  a    ^;)y)  ucx,o)  u.cx;o)  ax  =  o 

where     U»  (  X.o)  is  the  homogeneous  solution  to  Eq.    (32).    When 
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the  value  of       OCoa^ays.       >  from  Eq.   (29),  is  employed  the  function 


Ll    CXjO)        is  determined  to  be 


•      H  *     H  / 

U  Cx,o)     =      u «.     X 


The  value  of      (X       is  defined  by  substituting  the  functions        ^Cx)    , 

(fi(x)         (from  Eqs.    (13)  and  (18)),      LX(X,0)     ,  and      (JL     (X,o) 
into  the  solvability  condition  and  the  result  is 

Oi      -        ~     CX     .  (33) 


This  result  and  Eq.   (29)  determine  the  linear  approximation  to 
0^(0     ,  Eq.    (26).     Thus, 

a 


oCU)   -  (PA)    (/n'-zin2)    -    £ 


It  now  becomes  evident  that  the  question  posed  concerning  the  stability 
of  the  fundamental  mode  (  /V\—  X    )  for  sufficiently  small      €;       is 
answered  by  the    linear  approximation  for       q(  CO     •     Since  values  of 
/W\  >  1      give  very  large,  negative  values  of      C\  CO      and  thereby 
assure  the  stability  of  these  modes,  the  concern  is  centered  about 
/>Vi    "=  /y\    "=■  X    •    For  this  case 

ti  (O    =     -  €cl 

and  the  stability  of  the  fundamental  mode  will  depend  upon  the  sign 
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of  the  linear  term  of      Ca  (  €. )      .     The  various  cases  to  be  analyzed  are 
the  following: 

a)  negative  reactivity 

0((£)         <     0  for  positive        € 

>      O  ^°r  ne9ative       £ 

These  equilibrium  states  of  flux  are  depicted  in  Figure  2  and  the  corres- 
ponding eigenvalues  are  shown  by  the  solid  curve  in  Figure  3.    A  positive 
value  of      £        describes  an  equilibrium  state  of  flux  which  is  larger 
than  the  initial  flux.     This  state,  Eq .   (21),  exists  for  an  eigenvalue, 
A  (O  """    CX.      /  greater  than  the  initial  material  buckling  of  the  unper- 
turbed reactor,         A  o        »  an<^  ^  ^s  observed  to  be  stable.    A  negative 
value  of     £        indicates  an  unstable  equilibrium  state  in  the  reactor. 
However,  Eq .    (21)  indicates  that  the  eigenvalues  must  be  greater  than 
the  value  of     Ao     f°r  any  value  of     £       .     Therefore,  a  situation  in 
the  reactor  which  permits  a  decrease  in  flux  (i.e.  ,  negative     £.      )  with 
an  associated  increase  of  the  reactor  buckling  is  physically  unacceptable 

b)  positive  reactivity  (evaluated  by  changing  the  sign  of  all  terms  con- 
taining    01    ) 

O^  (  £  )  >         O  for  Positive    £ 

<         Q)  for  negative    £ 
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These  states  of  flux  are  depicted  in  Figure  2  and  eigenvalues  shown  by 
the  dashed  line  in  Figure  3  are  analyzed  in  a  similar  manner.    The  stable 
state  of  flux  is  that  described  by  a  negative  value  of     £       .    From  Eq.    (21), 
noting  that  the  sign  of   OL  must  be  changed  to  accommodate  positive  re- 
activity, the  fundamental  eigenvalue  is  less  than  that  of  the  initial 
material  buckling,      A©    .    Therefore,  a  decrease  in  flux  (i.e.,  negative 
t       )  is  compatible  with 'a  decrease  in  the  eigenvalues.    However,  for 
the  case  of  positive      £       the  equilibrium  state  is  unstable.    Again  it 
is  noted  that  the  increase  of  flux  associated  with  a  decrease  in  the 
material  buckling  is  uncompatible  and  this  situation  is  physically  un- 
realistic . 

The  solution  to  Eq .   (32)  can  be  completed  by  first  substituting 
Eq.    (33)  into  Eq.    (32).    The  right  hand  side  of  this  equation  becomes: 

RHS    -     [et(o)    +   a      «ft7x)J  LU'X/O)   • 

However,  from  Eqs .   (13),   (IS),  and  (30);         ^i*)    =     A/n    4/n  , 
(PCX)    =   /Iai  'Yin  j  and      UCXfi)  ~  IX^^T/m       respectively.    Thus, 

RHS   =     (-a  +   au>)  U^Y^ 

=    o. 

Equation  (32)  correspondingly  reduces  to  the  form 
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S7l  \lc%o)  +  ( Aom-  a  -  d0^)  uex,o)  ^  o  i*  & 


Thus,  the  solution  to         LA(X,0)       is 

Ut*,0)     =      UCX.O)    -      TX^^h     .  (34) 

For  larger  arguments  of     fc.       the  third  term  of  the  expansion  of 
OC  (  £  )      ,  Eq .   (2  6),  may  become  important.    This  term  may  be  evalu- 
ated by  differentiating  Eq.    (31)  with  respect  to     £        and  then  taking  the 
limit  as      £       tends  to  zero.    The  differentiation  yields 

y*  u(k,o  +  (  ^(t)  -  a  -  <^(0 


When  the  limit  as  £  tends  to  zero  is  taken  and  the  various  kno 
functions  are  substituted  into  the  resulting  equation,  the  following 
becomes  evident. 


wn 


36 


[it-    a(2-  fc)]  U^ 


(3  5) 


/v\  IW     5 


U(^O)      =      O  oyv     SO. 

The  solvability  condition  requires  that 

/  [&  -  a  (2  -/u)]  U^^utvad*  =0 

..    H 
where      UL  (OC^O)        is  the  homogeneous  solution  to  the  above  equation. 

»  H  „   H  .. .   H     ^ 

The  nontrivial  solution  for     U»  C%0)     is  LX  ("X,0)  =   v-*/fct       //»«*.• 


Thus, 


and  since  the  integra 


i^^d^l. 


/  ^  Cm   \ 

The  expansion  of     CX  (£)         ,  Eq.    (2  6),  may  be  written  in  the  more 
complete  form  by  the  substitutions  of  Eqs.    (29),    (33),  and   (36).     Thus 

+  K<*  (2-  §;)+••-.  (3, 

For  larger  arguments  of     fc        the  stability  of  the  system  is  contingent 
upon  the  negativity  of      0(    [c  )      t  defined  by  Eq.    (3  7)  .    For  the 
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fundamental  mode  (  /Y\  -  1    )  Eq.    (37)  reduces  to 

The  value  of   (YV\      of  particular  interest  is  that  given  by  fW\  -  1    since 
values  of  /YY\    >  1      provide  a  large,  negative  first  term  in  Eq.   (38). 
Thus,  the  stability  criterion  is 


or 


-£a(i-£+iec/A,)<o 

By  inspection  of  Eq.    (39)  the  quantity  within  the  parenthesis  may  be 
seen  to  always  be  positive  for  small  arguments  of      £       .    Thus,  the 
stability  is  governed  by  the  coefficient  ( -*  c  d    )  as  has  previously 
been  considered.    However,  for  cases  of  larger  arguments  of      £       the 
inequality  shown  above  may  be  changed.    Thus,  to  the  order  of        fc: 
the  condition  for  stability  is 
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1  -  i  7a 


(40) 


If  this  were  violated,  such  large  values  of     £       would  render  the  re- 
actor  unstable.    Thus,  the  undetermined  peak  amplitude  of      ^P  (  OO      / 
C,  ,  as  compared  to  the  peak  amplitude  of  the  unperturbed  flux,  A,  , 
is  important  when  considering  the  stability  regime. 
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C.         EQUILIBRIUM  STATES  FOR  A  FINITE  CYLINDRICAL  REACTOR 

The  finite  cylinder  represents  a  realistic  geometric  configuration 
for  a  reactor.    The  physical  and  mathematical  considerations  which  have 
been  deduced  for  the  slab  reactor  are  directly  applicable  to  the  cylindrical 
reactor.    In  this  section,  the  equilibrium  states  of  flux,  associated 
eigenvalues,  and  the  stability  of  these  states  which  have  undergone 
feedback  perturbation  from  some  initial  state  will  be  studied  for  cylind- 
rical geometry  with  azimuthal  symmetry.    The  flux  distribution  must 
depend  upon  two  spacial  coordinates  identified  by     \         ,  the  radial 
coordinate,  and      2?      ,  the  longitudinal  coordinate. 

The  equilibrium  solution  for  the  reactor  is  developed  from  the  govern- 
ing  equations  under  the  conditions  that  /     u*-  . 

The  governing  equations  reduce  to  the  form  given  by  Eq .    (10)  which  is 
rewritten  below, 

r-72  m,  /\  /)       C&*j2lI)  \    «  (41)    . 

t-Pir.-Ei*)    «    O  On     £}>, 

The  feedback-perturbed  solution  to  Eq .    (41)  is  evaluated  by  a 
Taylor  series  expansion  about  the  initial  state.     If  the  amount  of  per- 
turbation is       £     ,  which  will  subsequently  be  defined,  the  resultant 
equilibrium  state  is 
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Wr&t)    =    92(V*)  +   CWr,*)  +  it,    <P(rti)+>--   (42) 


I        *■" 


The  initial  state,       To^i  w        and       Ao      /is  defined  as  the  solution 
to 

^foCr,?)    +    (AG-a)    ftc*!*)    sO  in    D  <43> 

*/•(*"/*)  ■=  o  ok  £i>  - 

Equation  (43)  is  solved  by  assuming  a  separable  solution  for        toi^i) 
given  by 

<&(/;*)   -    R(r)  j?  <*;    .  (44) 

Substituting  Eq .    (44)  into  Eq.    (43)  yields  the  following: 

d    1  id/         dRCr)\ 

+  (Ao-a)RcD  get)    =   o 

which  may  be  separated  into  the  equations 

d     /      cAR(r)\  t  2  0 


d2  Zc*) 


-+      uj     2(^)    =o         u    j> 
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—      Z  *•  -N 

and  ^j  +        Co  "=         A0     -     OL     . 


The  solutions  to  these  equations  are  the  following: 

RCrl     =     A*        Jo  (/**%.)         » 

?c*>   =  A*     Cos  (^J/2  i) 


and 


>B-a  =  OVr.)*  +  rwi) 


where     /Am   are  the  roots  of  the  zeroth  order  Bessel  function.    The 
solution  to  Eq.    (43)  is  » 

=    Am  Jo  (/»«70  CosC^z  a) 

and        K-o,    =   (^-/f.)*  +  (Ta)1         (45) 

where  ^      _.      j^    3^  5^     .    .    . 

For  the  purpose  of  simplifying  the  subsequent  analysis  the  shape  function 
j/v\       is  defined  for  the  cylinder  as 

Thus, 


To/h    (^2:)       —       /i/y^       /^    . 


(4  7) 
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The  remaining  terms  of  the  expansion  of     ^f  (  If}  £,  £  J        and 
A  CO     '  9iven  by  Eq.   (42),  are  determined  by  differentiating  Eq.    (41) 
with  respect  to      £       and  then  taking  the  limit  as      £       tends  to  zero. 
These  steps  lead  to  the  equation 

^(r^)    ■=.    o  on    it,      (48) 

If       j   (1)  ^)  *s  to  have  a  nontrivial  solution,  the  solvability  con- 

dition (provided  in  Appendix  A)  may  be  envoked .    Thus, 

f  *  A  *         ^ 

/   -  "X  <£(*;*)  <P<o*>    2jy-  dr  d*  =  o      <49> 

where      (J>  /f  £j         is  the  homogeneous  solution  to  Eq.    (4  8).    This 


solution  has  a  similar  shape  function  as  Eq.    (43).    Thus, 

The  fundtions     j£  C*)  *t)  and     ^C^t)  ,Eqs.    (47) ,  and  (50) , 

respectively,  are  substituted  into  Eq .    (4  9);  hence, 

/   ~  ^  Am  'A,   ri  ^   2j  r   dr  d  2  -  O 


J) 

or 


-  >  Am  Bi  2j/  ^ 
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Since  the  integral 
I 


L 


im  rdr  dz    =    ^rl  J,ztym„) 


±       O  y 

the  solvability  condition  may  be  satisfied  only  if      A  =    O      •    This 
result  may  be  substituted  into  Eq.    (4  8)  and  the  value       y*  *'/^)       is 
the  solution  to 

Hence, 

</>(r,i)     =     <plrtu     =      Bm'i    .  <50*> 

The  parameter      £       is  chosen  to  represent  the  average  magnitude 

of  the  linear  term  in  the  expansion  of  Eq .    (42) ,  weighted  by  the  function 

■ — t — _7  . —   v       and  measured  in  units  of  maximum  initial  flux. 

TJ  r£     J,   (j*») 

Thus, 

\  J  r°  J> c/*^ 

Differentiating  Eq.    (51)  with  respect  to     t.      and  letting      £      approach 
zero  yields 
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A        =    D     f       ^    25"r  dr  d* 
=      B/n  • 

Therefore  the  linear  terms  of  the  expansions  of      j   (  ^  ^,  O     and 
>(€)     are 

(53) 
and  ^        ~        O  . 

The  third  terms  of  the  expansion  of  Vr  (lT  £,  6  )  and  ^  (£j 
are  determined  by  differentiating  Eq.  (41)  with  respect  to  £*.  twice 
and  letting       £       approach  zero 

S7z^P(nv    +  (Ao-a)  ^^2)   =  (54) 

a 


The  solvability  condition  requires  that 


/. 


( 


.  Z 


where       G?  (y^£)       is  the  homogeneous  solution  to  Eq.   (54).    The 

hi 

solution  for      (J)   ^  ^)        is  similar  to  that  obtained  from  Eq .    (43), 


thus 


44 


The  functions  y  (f/i)  and  7  (^  1 )  ,  determined  from  Eqs  . 
(47)  and  (53)  respectively,  and  (f)  (r  ?)  are  substituted  into  the 
solvability  condition.    Thus 

or      (a-x)  /,.f/J  j-drdz    -O. 


Therefore, 

A       "        CI       .  (55a) 

If  this  result  is  placed  into  Eq .    (54),  the  right  hand  side  will  identically 
be  equal  to  zero.     Thus  the  solution  for      7     (T  /  w     is  similar  to  that 
obtained  for  Eq.   (43)  and 


£<^>  =  cM  ii  . 


(55b) 
/VJ        'At 

The  various  results  for  the  terms  of  the  expansion  of  Eq .    (42)  are 
presented  as  follows: 


A(0-  Q-eVr.f+C^r  +  KA-*- 
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It  is  apparent  that  for  any  mode  of  flux  there  will  be  two  possible 
equilibrium  states  defined  by  the  negative  and  positive  values  of     £ 
For  each  of  these  states  a  unique  eigenvalue  is  present.    For  negative 
reactivity  (  Av^>)""GL)  must  be  greater  than     ^o        and  for  positive  re- 
activity the  opposite  must  be  apparent.     It  may  be  noted  that  the  eigen- 
values are  defined  by     X^'  """  01    •    Where  the  term       —   CL        denotes 
the  amount  of  change  of  the  material  buckling  to  accommodate  the  tem- 
perature feedback. 

D.         STABILITY  ANALYSIS  OF  THE  FINITE  CYLINDRICAL  REACTOR 

The  stability  criterion  for  a  cylindrical  reactor  is  developed  by 
considering  the  temporal  effect  of  a  deviation  on  the  equilibrium  states. 
The  governing  equation  has  the  form 


=  v'?W>  +(HVg)& 


a'c 

fCf,*,o)    =  <&  ens)   . 

The  temporal  nature  of  the  flux  is  investigated  by  imposing  an 
initial  deviation,    77        ,  on  the  equilibrium  states  defined  by  Eq.    (5  6) 
and  then  observing  the  growth  of  the  deviation.    An  unbounded  increase 
will  indicate  instability;  neutral  stability  implies  that  the  deviation 
tends  to  a  new  equilibrium  state;  and  asymptotic  stability  is  assured 
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by  the  decay  of  the  deviation  to  zero.  The  deviation  is  described  by  a 
Taylor  series  expansion  about  the  equilibrium  state  t  ^i  ^>^)  *  de- 
fined by  Eq.    (56): 

t-PCr^^.t)      =     ^<Xs,*J     +    y[  %  07VA>  (58) 

such  that 

'Vj— >o 

Equation  (57)  is  differentiated  with  respect  to     ^       and  the  limit  as 
"V?       tends  to  zero  is  taken.    The  result  of  this  process  is 

If  a  solution  for    y^j^^  ^j^i^'       is  chosen  to  be 

then  the  stability  of  the  system  is  determined  by  the  sign  of     (X(  w 
A  form  for    OC(£j       is  given  by  the  series  expansion 

o((0    =    oL0    +  £  oi    -t-x^^4*'-'        (59) 

The  value  of    OL  (c)    is  determined  by  the  perturbation  technique 
previously  used.     The  function  given  by  Eq.    (58b)  is  substituted  into 
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Eq.    (58a)  and  the  resulting  relation  is  ,„. 

V2u  +  (>(0  -  a  -  rite)-  a  j!^  t^J)^=° 

UL(r;  ^e)    =  O  on   a  , 

The  value  of     CX0      may  be  determined  by  taking  the  limit  of  Eq .   (60) 
as      £"       tends  to  zero.    Thus, 

z  /  \  (61) 

7   UOr;*)   +   (  Ae-  a    -  ola)  Ucr,  i)   =   O        in    £ 

A  nontrivial  solution  to  this  equation  is  obtained  only  if  (see  Appendix 
C  for  proof  of  this) 

XoW  -  a.  «   X.w  -  o.-  oCa   . 

Thus    OC  ©      is  given  by 

o(OMfm  =  (Ao*-  a)  -( "\om-  <x)  <52> 

=  ;k  (/u£-/i£)  +(%)  (/nl-/m'). 

Some  remarks  may  be  directed  toward  the  stability  of  the  system  via 
Eq.    (62).    As  was  apparent  for  the  slab  reactor  the  fundamental  mode 
(  f/\  -z.  1   )  may  be  evaluated  as  follows: 

Oioyrrr,        =      O  ff*     =   1  (62a) 

<      -  2  TJ  Av>     >    1     . 
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Thus  the  higher  modes  of  deviation  (/W\>1  )  which  are  imposed  upon 

the  equilibrium  flux  will  vanish  with  time.     However,  the  value  of    /yv\  ^r.  ^ 

provides  0\o\m\^:'  ^  anc*  t*ie  system  may  ^e  thought  of  as  neutrally 
stable.  Thus  it  is  necessary  to  inquire  about  the  effect  of  OL  upon 
the  system  for  the  value  of   fYY\  —  1  . 

The  higher  modes  of  the  reactor  (  f/\  >  j[   )  are  analyzed  as  follows: 

oi(^)     <    -2tj-  /w>    >  /n 

=  O  /vr\     =  m*  (62b) 


> 


2.  tj  on  <  /n.  . 


Since  the  condition  provided  by    /W\    ^  /V\       may  be  realized  for  modes 
of    fY\   ">  x.      ,  values  of    C\(£>^  O  are  present  and  the  higher  modes 
are  unstable.    However,  the  reactor  is  assumed  to  be  operating  at  the 
fundamental  mode  before  the  feedback  is  introduced.     Equation  (62)  for 
0\  q(y\  m\      may  be  used  in  Eq .   (61)  to  provide  a  result  similar  to  Eq .   (43)  . 
Thus, 

lACr,*)     =      UL/m     im*  (G3) 

(See  Appendix  C  for  the  proof  of  this.) 

The  linear  term  of  the  expansion  of    0\  (£)   (Eq .    (59))  is  evaluated 
by  differentiating  Eq .    (60)  with  respect  to      £       and  taking  the  limit  as 
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C      approaches  zero.    The  result  is 

y\x(r;*)   +  ("Xo^  -  a  -  o(0r,J   u(ir;i)      (64) 


By  noting  that     /\       has  previously  been  determined  to  be  zero  the 
solvability  condition  requires  that 

f     f   \  tyftl 

J  r 


[c<  +  a  ^t^jJ  uku  u^*)    =  o 


where        (.X   (r(^)      is  the  homogeneous  solution  to  Eq.    (64).     By  the 

same  reasoning  used  for  the  solution  of       LL  (^  :k)    the  solution  for 
,     N 

UL  tr,2)    is 


U-Cf,S)      =       LC^      %^ 


When  this  result  and  those  provided  by  Eqs  .    (4  7),    (53),  and  (63)  are 
used  to  evaluate  Eq.    (64b)  a  value  for     (X       is  derived.    Thus, 

or 

(ot  +  o.)]    i*m  r  dr  dz    =  O  . 
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Therefore, 

© 

oi  =    -  a .  (65) 

This  result  may  be  substituted  into  the  expansion  for      CX  \€)       and 
further  analysis  of  the  stability  is  made.    Thus, 

0(  (£)       =       oComrm      -     £<X     .  (66) 

The  question  pertaining  to  the  stability  of  the  fundamental  mode  ( /v\  =  J_) 
is  resolved  for  small      £.      .    For  the  case  where    /yy\  =  1     the  value  of 
(X  o  1 1  is  zero  and 

oite)     =     -So.. 

The  analysis  of  the  various  cases  of  negative  and  positive  feedback  is 
similar  to  that  developed  for  the  slab  reactor.    The  conditions  of  stability 
for  the  fundamental  mode  are  briefly  stated  as  follows: 

a)  The  stable  equilibrium  state  for  a  reactor  with  negative  reactivity 

has  a  positive  value  of     c       •    The  other  equilibrium  state  is  unrealistic. 

b)  The  stable  equilibrium  state  for  a  reactor  with  positive  reactivity  has 

a  negative  value  of      t      .    The  other  equilibrium  state  is  unrealistic. 

If  the  value  of    OC      is  substituted  into  Eq.    (64),  the  right  hand  side 

will  become  zero.    Therefore  the  solution  for        Lk.  OG  w       is  the 

homogeneous  solution  to  Eq .   (64)  and 

UL   (  *",  l)       =r         IX  rrv\      J  /y^    -  (66a) 
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For  larger  values  of     c       the  quadratic  term  of      (A  (w      may  be 
necessary.    To  obtain  the  desired  results  Eq.    (60)  is  differentiated  twice 
with  respect  to      £.       anc*  the  limit  as      £       tends  to  zero  is  evaluated. 
If  all  of  the  known  functions  are  employed  the  result  is 

V      U«;*)       +   (Xo/n    "    0.     -    Oiommv)    iliXi)      =(67) 

(<£-  a  (2-C"/aJ)  TJU  "k*       in  * 

u.  c^d    =  o  ova     s^  . 

The  solvability  condition  requires  that 

|[o(  -  a(2-^)]u^ALuHCY;^  2jrdrd^  =0 

,,    H 

where     IX  £/*,£)  is  the  homogeneous  solution  to  Eq.    (67).    Therefore 

"  tf   _TTH    Ai/ 

UL      "*  LA^    T/w     •     It  is  obvious  that  the  solvability  criterion  is  met 

only  if 

oL  -    a(2-  fe)     =  O 


This  result  is  substituted  into  the  expansion  for     (X  C&  )      ,  Eq.    (59), 
and  the  following  becomes  apparent, 


(68) 
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(69) 


-€a  +{f&(2-C*A„) 


Since  the  fundamental  mode  of  operation  is  generally  of  greatest  interest, 
Eq.    (69)  is  rewritten  for      fA  ~  1     . 

rf<0    «  Or;)  (p^-^m)    +("%)  (l-/m)<70> 

-  £a[l  -  l€  (2-c/Al)]  +•••. 

For    /WJ   >  1      the  first  two  terms  dominate  and      0\  (w  *^  O    •    For 

/Vv\  =  1    the  stability  may  be  determined  by  evaluating 

oC(«)  =   -€a[l  -  i€  ^2-VaJ]- 

These  results  are  similar  to  those  obtained  for  the  slab  reactor  and 
stability  is  granted  under  the  conditions  established  for  Eq.    (37). 
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III.     ANALYSIS  OF  THE  REACTOR  WITH  A 
TEMPERATURE  COEFFICIENT 

**-=.-  a  T"3/2     (T"3/2  MODEL) 


The  solution  of  the  nonlinear  equations  (3),   (4),  and  (8)  governing 
the  reactor  kinetics  is  obtained  by  the  same  perturbation  technique  used 
to  analyze  the  "T       model."    This  method  investigates  the  effect  of  a 
small  perturbation  introduced  by  the  feedback  effect.    The  equilibrium 
state  that  ensues  from  the  perturbation  is  defined  by  a  series  expansion 
about  the  initial  unperturbed  state  of  the  reactor.    The  perturbation 
technique  offers  multiple  solutions  and  unique  eigenvalues  for  each  of 
the  equilibrium  states  (see  Figure  2).    The  presence  of  more  than  one 
equilibrium  state,  known  as  the  bifurcation  phenomenon,  requires  an 
investigation  of  the  stability  of  each  of  these  solutions.    The  procedure 
of  analyzing  the  stability  entails  imposing  an  initial  deviation  on  the 
equilibrium  flux  and  then  examining  the  resulting  temporal  behavior  of 
the  deviation. 

A.         EQUILIBRIUM  STATES  FOR  THE  SLAB  REACTOR 

-3/2 
The  first  consideration  to  be  made  in  the  analysis  of  the  "T 

model"  is  that  of  determining  the  possible  equilibrium  states  of  the 

reactor  after  a  feedback  perturbation  has  been  imposed.    The  solution 

to  the  governing  Eqs.    (3),    (4),  and  (8)  is  desired  for  the  situation  where 


th 


e  perturbed  flux  has  achieved  the  equilibrium  state  (i.e.  ,       7\  ^     ~®) 
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The  governing  equations  reduce  to  the  form 


and 


The  equation  for  the  multiplication  factor  is  easily  shown  to  be 

The  function     y^  (X)       and  the  constant     ^0       are  considered  the 
initial  flux  distribution  and  multiplication  factor  before  the  feedback, 
effect  becomes  significant.    The  result  is  substituted  into    Eq.    (3), 
yielding  thereby: 


where 


A    =    L  -  1. 


After  the  perturbation  due  to  feedback  the  new  equilibrium  state 
is  described  by  a  Taylor  series  expansion  about  the  initial  flux  distri- 
bution,      /o  £  **  )      i  and  the  initial  eigenvalues,     /\  0      .If  the  amount 
of  perturbation  is  measured  by     ^      (the  perturbation  parameter  which 
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will  subsequently  be  defined),  the  resultant  reactor  state  is 

'X(t-)    =    Ao    +   €  A    +  {£  A  +  '•  '  . 

For  small  perturbations  of  the  system  the  values  of    j(x,£)    and     ^  (O 
are  described  by  truncating  the  series  after  the  third  term.    The  un- 
perturbed states  for  the  reactor  is  defined  as  the  solution  to  the  problem 

V^CX)         +       )\Q    <&(X)        =     O  in   D  (72) 

<#>(*)      =       O  ow      ciD  , 

It  is  readily  seen  that 

/n     •=    l ,  J,  s,  •  •  •   .  -/, 

So  as  to  simplify  the  subsequent  analysis  the  shape  function     T/h    is 
defined  for  the  slab  reactor  as 

4m   -   Cos  ("Vz)  x  .  (73) 


Thus, 


<&  Cx)     =     Am   "/%  .  (74) 


The  next  term  of  the  expansion  is  determined  by  differentiating 
Eq.  (71)  with  respect  to  £:  and  taking  the  limit  as  £  tends  to 
zero.    This  differentiation  yields 
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r  'Xt.  i    (75) 

=   ~  (  /\&)  -    I    V^ft)  fivj)  *PC\ti     'in  J) 

As      £       approaches  zero,    T^C'Kftz)        approaches      To^J       ,     A  (^^ 
approaches    A©         /  an^ 

In    i) 

In  order  to  obtain  a  solution  to  Eq .   (7  6)  the  solvability  condition  (given 

in  Appendix  A)  requires  that  the 

r  u 

/    (-^Qcx))    fit*)     dx     =  O  (77) 

where     r  Cx)      is  the  homogeneous  solution  to  Eq.    (76)  given  below 

z    .    H  r  0.   ^-/i  -|       •   //  (78) 

4^    <X)  «       O  OKI  <fJ>. 

A  solution  to  Eq .  (7  8)  is  offered  by  a  modification  to  the  R.  Courant 
and  D.  Hilbert  perturbation  method  (provided  in  Appendix  B) .  Their 
approximated  solution  is  given  by 

•   n  a  /a\l  (79) 
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(ca/  \L 

Since  the  contribution  of  \  /  2.  /  is  very  small  only  the  first  two  terms 
of  the  series  will  be  retained.  The  substitution  of  Eq.  (79)  into  Eq .  (78) 
leads  to  the  requirement  that    Ll^(X)     be  the  solution  to 


Um(X)       -    O 


ovx       SJ> 


the  solution  to     A/m()()     can  be  obtained  as 


oo 


aj^(x)  =  y 


CI  /n 


A  om    ~~   A  oa  J 


where 


d 


rv\ 


I  «2 


;    Us*  U;    dn< 


■=  <t 


The  orthonormal  solution  for    U./^X)      is    U./i^Xj 
Ae^        ~      (        '-yz  )  •     The  constant      Ca /y\l 


av       and 


is 


determined  as 


d 


MJ 


=  /7C 


/  ^  ^   dx 


thus 


oo 


d 


V  (if  — ^j oh. 

These  results  may  be  combined  into  the  expansion  given  by  Eq.   (77)  to 


yield  the  following  result, 


oo 


d 


(80) 
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For  the  convenience  of  the  subsequent  analysis  this  function  will  be 
referred  to  as  the     ]_      function  given  by 


r„  =  <  +^£  -^h-^ 


AV 


J  ±  /K. 


J 


-f    .    .     . 


(si; 


A  question  may  arise  concerning  the  convergence  of  the  series 
used  in  Eq.    (81) .     If  the  Weierstrass  M  test  for  uniform  convergence  is 
employed,  a  series  which  is  a  continuous  function  of     A       in  the  domain 


D  is  convergent  under  the  condition 


i 


where 


Em 


is  a  known  convergent  series  of  positive  con- 


stants [6]  .     If  the  convergence  of 

.  d/i 


£<*>  -  '%??  i«> 


is  doubted  this  test  may  be  used.    First,  the  maximum  value  of 


is  determined  as  follows: 


+1  u 


-TV"'"  -y 


i/vi  -1  -' 


dsr 


v^/-,    dx 


6     ^ 
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Since  the  function       J  j 

2 


/7?„  (j'-m) 


Y 4     -~  I      throughout  the  domain, 


d  ma 


J 


If  the  convergence  of  the  left  hand  side  of  the  inequality  is  provided, 
then    |  *  ix)       is  convergent  by  virtue  of  the  Weierstrass  M  test.    The 
integral  test  provides  a  convenient  method  to  establish  this  convergence 


Hence,  if  the  integral 
CO  -j 


/ 

J J 


z  .  i\    d  A 


JK(f-ri)    UJ 


exists,  its  convergence  and  that  of      - — 7 J    •  i        ijCXJ  is 


determined . 


Thus, 


J3  &Z  (j1-^)  dj    =  jfc(^)Co+h^ 


CO 


2-./-L' 


l  -  CoHJu) 


where      J        is  any  possible  value  of      j        .     Since  the  integral  exists, 
the  series  converges. 

The  values  of      %  (*)       and    ty  (  Xj         from  Eqs  .    (74)  and  (80)  , 
respectively,  are  substituted  into  Eq .    {77)  and      /\       is  evaluated  as 
shown . 
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/- */u ± sn  \±* 2#y  ■£&? 4;+" ]dx-o. 

J»  J   few-       J        J 

Because  of  the  orthogonality  of  the       7^        function  the  above  equation 
reduces  to  the  following  form, 

M/UB^  ^    dx   =0. 

This  equation  may  be  satisfied  by  the  condition 

■\     =     O.  (82) 

The  value  of     /|        is  substituted  into  Eq.   (76)  and  the  homogeneous 
equation  that  results  is  similar  to  that  defined  by      ^r  ^  <*•)      •    Thus, 

¥>CX)      -      B„      T*.  (83, 

The  value  of  the  perturbation  parameter,     £       /  is  chosen  to  represent 
the  average  magnitude  of  the  linear  term  in  the  expansion  for  Eq .    (11) 
weighted  by  the      X/y,     function  and  measured  in  units  of  the  maximum 


initial  flux.     Thus 


Vd^xilVw-*."*]**     (84) 


If  Eq.    (84)  is  differentiated  with  respect  to     £       ,  the  result  is 
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By  letting      £"       approach  zero  and  employing  the  solution  for      ^P (x) 
(i.e.,  tp(X)  -  U/y\  L  )  the  above  equation  reduces  to 


t  -I,  t:  ** 


or,  from  equation  (81) 


A 


M 


& 


M 


/J 


x  +  j£  &r  * + 


'J 


] 


dtf 


The  quantity  within  the  parenthesis  is  expanded  and 


A/A        __ 

B 


Al 


d%    + 


d  /v\ 


/y\    - 


•i       h»'idx 


+  a.-/(§i)' 


Do 


,.j=lj3,S,- 


T  J~  'J 


d*  +  0«?). 


The  second  term  on  the  right  hand  side  is  identically  zero  because  of 
the  orthogonality  of  the       iM       function.    When  the  integration  of  the 
first  and  third  terms  is  evaluated,  i.e .  ,      I -^      l/y^  Q  X  ~*  /  the  result 


is 


Am 


=    1  +  a 


-i2 


I 


+ 


0<a') 


^  At. 


_  4. 

Since  the  order  of    CL        is  generally  very  small,  the  contribution  of 
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the  higher  order  terms  are  neglected.    Therefore, 

m 

Equation  (83)  that  defines      ^r  ( "X  )        may  be  rewritten  as 

<?t*)    =    A*    ?„  <85> 

The  third  terms  in  the  expansions  for    yl (x  &)    and         \  (  6) 
given  in  Eq.    (11),  are  determined  by  differentiating  Eq.    (75)  with  respect 
to      £,       and  then  taking  the  limit  as      £       approaches  zero.    The  dif- 
ferentiation yields 

+  ft(0     -      2      ^tx/fej    <£<x,6>)        V<X,6)      = 


As     £        approaches  zero  and  Eq.    (82)  for     /\CO)     is  used  the  above 
equation  simplifies  to  the  form 
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/ 


V*#(x>    +    (*0    -     7  ^0"cxL))  9<x)      =     <86> 


a. 


>     foCX)       +    %    4*<.x)    <P  CX)         ini 


The  solvability  condition  for  Eq.   (86)  is 

where      (£>(x)       <  the  homogeneous  solution  to  Eq.    (86), may  be  seen 
to  have  precisely  the  same  shape  function  as  the  solution  to  Eq.    (78). 

ThuS'  H 


If  the  values  of     9^  (x)        ,     *P(X)        ,  and    *P (  X)       from  Eqs  .   (74), 
(85),  and  (88),  respectively,  are  substituted  into  the  solvability  con- 
dition of  Eq.    (87),  the  following  becomes  evident, 

The  integrals  are  simplified  using  the  definition  of       J_.  a\      given  in 
Eq.    (81).     Thus 
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•  ft 


CO 


2ar     ^/ 


u 


dx 


dx 


-    1. 


because  of  the  orthogonality  of  the        //>?        function  and 


* 


/j7o^x=/k*+f, 


^M.i_  4 


3—  _D_ 


CO 


n3 


dx 


If  the  contribution  of  any  term  having  a  coefficient  of     CL     or  larger  is 
considered  negligible,  an  approximate  value  of     /\       may  be  evaluated. 


Hence , 


H 


a/uc:  = 


Stiill^d* 


(89) 


The  definite  integral  may  be  evaluated  for  the  fundamental  mode  of 
reactor  operation  (i.e.  ,   /Vl=  1    )  and  the  value  of     /\        for    /V\  -  i    is 


a 


X     a     ^   (-279)  . 


(90) 


The  equilibrium  solution  for     j(Xj6r)     and      /^  (c/         is  obtained  fron 
the  results  of  Eqs.    (72),    (74),    (82),    (85),  and   (90).     The  expansion  for 
tfCX^)     and       ))  (£)        ,  equation  (11)  is 


^<%«  =  AmI^*  £  71  +  o  a1))  . 


or 


65 


where        i/v\        and       X/Vi       are  defined  ^y  Eqs.    (73)  and  (81),  respec- 
tively, and  for  the  fundamental  mode 


=  (%)*    +     -03^9^    f   ---.(92» 

The  magnitude  of  A.   is  determined  by  the  "power  condition" 

p  -    o  I    cp(X)  doc 

where      r        is  the  reactor  power  level  and      Q       is  a  conversion  factor. 
Equation  (91)  implies  that  there  are  two  equilibrium  states  for  each 
mode  of  operation  corresponding  to  the  negative  and  positive  values  of 

C      •    For  the  fundamental  mode  this  condition  may  be  depicted  as 
shown  in  Figure  2.    For  reasonable  values  of      £       and  ( (X.  )  the  con- 
tribution of      _L  /v\      is  essentially  that  given  by        7  a\  The  impact 
of  the  series  of  cosine  terms  within  the      J_  /rt      function  may  be  seen  to 
be  negligible  if  the  peak  amplitude         C      ri  13     MAX  of  the  largest 
contributor  is  considered. 

2_a       d i3        ^l   I 

.WHH 

The  value  of      0*  ,->    —  77T7  was  calculateci  anc^  f°r  a  very 

^  '     A  Q- 

conservative  value  of         c  ~  •  1       an^     IK*       ~    .  OO  J. 
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c  H  I3M/9X  -J-  ir2-  i  8       ^ 

«     2.91    x   io"6 

«  1  . 

For  each  of  the  equilibrium  states  there  are  unique  eigenvalues  as  shown 
in  Figure  3.    For  a  reactor  with  a  negative  coefficient  of  reactivity  the 
fundamental  eigenvalue  is  drawn  as  a  parabola  centered  at    (     -y  2. ) 
which  opens  to  the  right.    Thus  the  equilibrium  eigenvalues  are  always 
larger  than  those  of  the  initial  state.    As  has  already  been  deduced  for 
the  "T       model"  any  decrease  in  the  neutron  flux  (i.e.,  negative     c.       ) 
cannot  be  compatible  with  an  increase  in  the  material  buckling.    There- 
fore negative  values  of     £        for  negative  reactivity  are  unrealistic. 
For  the  case  of  positive  reactivity  the  eigenvalues  react  in  an  opposite 
way  as  depicted  by  the  dashed  line.     Values  of      £      which  are  negative 
represent  possible  states  of  flux,  but  positive  values  of      £       are  un- 
realistic. 

When  the  eigenvalues  of  the  "T       model"  are  compared  with  those 

-3/2  -1 

of  the  "T  model"  on  Figure  3  it  is  noted  that  the  "T       model"  is 

affected  more  by  equal  perturbations  on  the  system.    This  result  is  con- 
sistent with  that  predicted  by  A.  Reichel  [3]  as  can  be  shown  on  Figure  1 
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There  it  is  observed  that  the  line  for  the  "T       model"  is  in  very  close 
proximity  to  Reichel's  line  of  maximum  Doppler  effect.     Hence,  one  may 

conclude  that  the  Doppler  effect  has  a  more  pronounced  impact  on  the 

-1  -3/2 

"T       model"  reactor  than  on  the  "T  model"  reactor. 


B.         STABILITY  ANALYSIS  OF  THE  SLAB  REACTOR 

The  stability  criterion  for  the  equilibrium  states  of  a  slab  reactor 
with  negligible  concentrations  of  delayed  neutrons  is  developed  as 
follows.     The  governing  equations,  Eqs  .   (3),    (4),  and  (8)  are  combined 
into  the  equation 


■b* 


=  V  Ww  +[X-  <x(<P~£v 


The  stability  of  an  equilibrium  state  is  investigated  by  prescribing 
an  initial  deviation,    *Vf       ,  in  the  equilibrium  flux  and  observing  the 
growth  of  this  deviation.     If  the  effect  of  the  deviation  increases  the 
flux  with  time,  the  system  is  said  to  be  unstable.     If  the  effect  is 
reversed  the  system  is  stable. 

The  time  dependent  flux  is  described  by  a  Taylor  series  expansion 
about  the  equilibrium  state,      /  t  ^  6)     ,  defined  by  Eq .    (91): 
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(P(Xjt,i:,yi)   =  ^(^k)   +  y{  41  (x,*,*)  (94) 

such  that        {j'Zof  <*>*.*<*)    *   <ftxy*). 

Equation  (93)  is  differentiated  with  respect  to    "V7       and  the  result  is 


-  <?*<£,  *  [>«  ♦  a (*¥*  <e&)]  4>   k> 


where 


^  «   o 


Ov\         <£j> 


If  the  limit  as     "W       tends  to  zero  is  taken,  the  equation  reduces  to 

-h        .*-y%i  .*   (95) 

in  J> 


l|.  »^v,+[aco  -  h;V  a #£]<{> 


^    =  o 

The  solution  for      iC.V'X^  ^  )         is  chosen  to  be 


ox 


ij). 


^Oc^tj    =     uck^j  e#p  [u (&)?]. 


(96) 


Hence  the  stability  of  the  system  may  be  analyzed  by  inquiring  whether 

OC  (£)      is  negative  or  positive.    The  value  of      cLL&.)     is  described 

by  a  Taylor  series  expansion  about  the  parameter      £      and  is  given  by 
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If  the  proposed  solution  for     <^L  (^T/  t)       (i.e.  ,  Eq.    (96)  is  substi- 
tuted into  Eq.    (95),  then 

vV(x;t)   -f  [^j  -o<(0  +  f  W*,o     (98) 

U*  C  x,  O         =  O  O*        £  J>    . 

Taking  the  limit  as      c       approaches  zero  results  in  the  expression: 

_-*  /-,  a,  //5  -}'*  \  (99) 

U(xyo)      =•    O  on     «£D   . 

A  solution  to  this  equation  is  offered  by  a  modification  to  the  R.  Courant 
and  D.  Hilbert  perturbation  method  (provided  in  Appendix  B) .     Their 
approximated  solution  is  given  by  the  series 

—  5*     -=.,  ,    /5r\Z    — -  (100) 

/a/  \L 

Since  the  contribution  of    V.   /2»  /       is  very  small,only  the  first  and 
second  terms  of  the  series  are  retained.    The  technique  provides  that 
(J-  (■  *)       is  the  solution  to 

V*ucx)     -t    (/\oon  -dj)   CLcx;   ^O      i^D 
U(  x)     =  o  ©*     <^£ 
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Since  the  value  of     A  om     nas  already  been  determined  in  Eq .   (72), 
the  only  nontrivial  solution  for    LA.  (  X)      is 


*)        -         W/wt     M 


where        f/yn        is  defined  by  Eq.    (73).     (Proof  of  this  solution  is  pro- 
vided in  Appendix  C)  .    The  eigenvalues  corresponding  to     LA ^i  *)  are 
Ao/yy\     ~      A  o/v\  ~*    <^\o      •     Therefore  the  first  term  for  the  series  for 
0(  (£)      is  defined  as 


Oiorv\/yn  Aorv\    ~     A 


onm  4 


The  linear  term  of  the  series  of  Eq .   (100)  is  given  by 


ft  /m 


where 


The  solution  for     LA.      may  readily  be  seen  to  be 

d, 
u=UM,IX»+  TV    TT-'T"  T/w  + 


where 


(101) 


This  solution  may  be  put  into  the  form  afforded  by  Eq .    (81) 
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If  Eq.  (101)  is  considered,  some  conclusions  as  to  the  stability 
of  the  reactor  for  very  small  €  can  be  made.  For  the  fundamental 
mode  of  the  reactor  (  /y)  -  J[     )  the  value  of     (Xoi/m     is 

CX  01  ivr\       —         A  o\       ~~     A  o  a**. 

-  (Wz)z  (l-v)  . 

The  possible  values  of    o(a\/rr-      are 

Olo\nn  "=        O  /ym.     =    1  (102a) 

^     -  2jj-  /m    —  3,  S,  7,  • •  . 

Thus  fcr  values  of       £        near  zero  the  first  term  of  the  series  for 
£X  (£/       indicates  that  the  flux  is  neutrally  stable  for  the  fundamental 
mode  of  the  deviation  (/VKI-  \  )  and  stable  for  the  higher  modes  of  the 
deviation  (/VH  >1  )  .    The  question  of  neutral  stability  obliges  an  in- 
vestigation into  the  remaining  terms  of   £\  (  O  . 

Some  conclusive  results  may  be  obtained  concerning  the  higher 
modes  of  reactor  operation.    The  possible  values  of    o(  (&)    for  the  case 
where    /V?   >  j.    are 

z=      o  />^7  =  /y\         (102b) 

2. 

^     2  tj  An  ^   /n  . 


a 


Since     o(  ((:)    may  be  negative,  the  higher  modes  are  unstable. 

The  linear  term  in  the  series  expansion  of      Os(t)    is  evaluated 
by  differentiating  Eq.    (98)  with  respect  to      €.       and  thus 

*  UCX,£)       =    O  . 
The  value  of       A(^)    has  already  been  determined  to  be  zero.     If  the 

limit  as     £        approaches  zero  is  taken,  then 

V      U(\C)      +    [Ao    "     2    **©         -    ^o/virtrj    U(XyO)    (103) 

/a    n~*A      •  .  V 


The  solvability  condition  requires  that 

f     (  ^CD~1/Z    ld  v\  •    «         ,  (1°4 

J      (    q    %<**)   Vex)      +    oCj    ucx)  LKCX)  dX      =0 

.     H 

where    La.        is  the  homogeneous  solution  to  Eq.    (103).     If  the  value 

°f     ^Oto/>n     is  substituted  into  the  homogeneous  equation,  the  result 
is  an  equation  similar  to  Eq .   (99),  thus 

^/m  LA   rm       L/yk  ' 

•     H 

When      LA  /vn         and      iX  An-         are  used  to  evaluate  the  solvability 

condition  the  value  of     0\      is  developed.    The  details  of  this  procedure 
are  developed  by  using  the  results  for       ji  (X  J       /      T    (  X)        >  anc^ 
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tC    (X)        obtained  from  Eqs  .    (74),    (85),  and  (102),  respectively 
and        U>    (y)        .     Substituting  these  functions  intoEq.   (104)  gives 


I  f?  C  i?l?M  +  <k)  ux,  ul  r„d<x~o 

or  ? 

Consider  the  coefficient  of     0(       in  view  of  the  definition  of        J.  , 


given  by  Eq.    (81) 


/rid*  -I  fa  *  f£&r  1& 

=    1  +  O  +    <9<ra') 
=    1. 

An  approximation  for  the  linear  term  of  Eq.    (105)  may  similarly  be 
obtained  to  be: 


+  o  m 

if  terms  of  the  second  power  of  (  CX.   )  and  larger  are  considered  negligible 
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These  approximations  are  substituted  into  Eq.    (105)  and 

(107) 


A  =  -i^/tVjd*. 


The  stability  of  the  system  is  determined  from  the  evaluation  of 

by  substituting  the  results  of  Eqs  .    (101)  and  (107)  into  the  expansion 

of  Eq.    (9  7).    Hence, 

otto  =(%)V-  mx)  -  4  €  kL'+f'K  *«  +""' 


The  additional  term  enables  the  further  discussion  of  the  stability 
of  the  fundamental  mode  of  reactor  operation.     Thus,  for       /Y\  —  1 


8) 


*  J       (109) 


For  all  values  of   /YV\  >1    the  first  term  becomes  dominant  and  forces 
the  value  of    0((£)     to  be  negative.    However,  the  choice  of     /YY[  —  J. 
defines  the  value  of       (X  (  w     # 


tffc>-  o-^Jk.ff*  d* 


This  integral  may  be  evaluated  and  the  result  is, 

diie.)  =    -€   pff  ('•279)  . 


This  result  provides  the  following  conclusions; 


(110) 


(111) 
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a)  negative  reactivity 

C>(    (0  <     O  for  positive      £ 

£>((£)  ^     O  for  negative      £. 

These  results  indicate  that  the  equilibrium  solution  described  by 
an  increase  in  the  flux  over  the  initial  flux  distribution  will  be  stable. 
On  the  other  hand,  the  other  equilibrium  solution,  which  was  determined 
to  be  unacceptable,  is  unstable. 

b)  positive  reactivity 

0(  (€)  >      O  for  Positive     €. 

q(    (£)  <     O  f°r  negative    £ 

The  equilibrium  state  described  by  a  negative  value  of       £ 
is  observed  to  be  stable.    The  other  equilibrium  state  which  was  con- 
cluded to  be  unrealistic  is  unstable. 

In  summary,  the  stability  of  a  nonlinear  reactor  was  studied  by 
solving  the  boundary  value  problem  and  then  considering  the  effect  of 
an  initial  deviation  on  these  solutions,  called  the  equilibrium  states. 
The  possible  equilibrium  states  are  those  described  by  Eq .    (91)  with  a 
unique  eigenvalue  for  each  state  given  by  Eq.    (92).     However,  the 
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stability  criterion  predicts  that  only  one  of  these  states  are  stable  and 
that  state  will  be  physically  observed  in  the  reactor.    For  a  negative 
temperature  coefficient  the  stable  state  is  one  described  by  a  positive 

£     .    For  the  case  of  positive  temperature  coefficient  the  negative 
value  of     £       represents  the  stable  condition.    Other  equilibrium  states 
are  unrealistic . 

C.         EQUILIBRIUM  STATES  FOR  A  FINITE  CYLINDRICAL  REACTOR 

The  analysis  of  the  equilibrium  states  of  the  finite  cylindrical 
reactor  can  be  done  in  a  manner  similar  to  that  done  for  the  slab  reactor. 
The  equilibrium  states  and  the  associated  eigenvalues  will  correspond 
to  those  developed  for  the  slab  reactor.    However  differences  will  be 
noted  because  of  the  differences  in  the  geometry.    The  flux  distribution 
depends  upon  a  radial  coordinate,     I  ,  which  varies  from  the  center 

line  to  the  outer  radius,       \0         /  and  a  longitudinal  coordinate,    i?      , 
which  varies  from  -1  to  +1. 

The  equilibrium  solution  is  developed  from  the  governing  equations 
which  are  simplified  by  requiring  /  Tit*    """    (—^   - 

The  governing  equations,  Eqs.    (3)  and  (8),  reduce  to  the  form, 
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The  perturbation  solution  is  evaluated  by  a  Taylor  series  expansion 
about  an  initial  reactor  state,     TqC^jH)     and         /\o       •     H     £       is  the 

perturbation  parameter,  which  will  subsequently  be  defined,  the  resultant 
equilibrium  state  is 

(pcwi   =  #cc»  +  e  ^<tt*>  +  z^ficm)  +>■  ■  an) 
)iai    «  *0   +  £  x   +  2t*  X  +  ••  • . 


The  initial  states,      tL^Z)    ,  of  the  reqctor  is  defined  as  the  solution 
to  the  system: 

The  function     To^^)       is  determined  in  a  manner  similar  to  that  already 
provided  for  Eq.    (43).    The  result  is 

<&(m>  =   Am  J0  (/*%-.)  Cos  C*¥z  z)  . 

ion         T  /t\       i 

%  -    JofcS/r.)    Cos(MJ/z  2) 

The  eigenvalues  are  also  determined  by 

where      yt^/w      is  the  n      root  of  the  zeroth  order  Bessel  function. 


If  the  shape  function         « /n       is  defined  for  the  cylinder  as 


(115) 
(116) 

(117) 
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The  linear  terms  of  the  expansion  of  Eq.    (113)  are  determined  by- 
differentiating  Eq.    (112)  with  respect  to     €  ■     and  then  taking  the  limit 
as      £       tends  to  zero.    The  result  of  these  operations  is 

S7Z</>(r;V    +  (A*  "    f  <&"£*))  $cr,*)  =    ("8) 

to 

In  order  that  a  solution  is  obtained  for  Eq.    (118)  the  solvability  con- 
dition (given  in  Appendix  A)  requires 

I   _^</>o(;v:2)  cpcr,*)    2rjr  dr  d2  =  O  (119) 

H 

where        ^~P  if  1)  *s  tne  homogeneous  solution  to  Eq.    (118).    The 

•    // 

equation  defining        y?  Cf,  2  )       /  Eq.    (120),  is  solved  by  a  modification 

to  the  R.  Courant  and  D.  Hilbert  perturbation  method  (provided  in 
Appendix  B)  .    Thus, 

.20) 


The  approximated  solution  is  given  by  the  series 


yV?aj    =    u^r^zj  +  f/v^C^a)  +  O^cf)  .  u«) 
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Only  the  first  two  terms  of  the  series  are  retained  since  the  contribution 

1 
of     Q-       is  small.    The  solution  for       U.^  ( *7  2: )      is  obtained  from  the 


equation 


V2uM(r^^     +    X0  UmOr,tL)     =   O  in  T> 


Since  the  solution  for  this  equation  was  already  determined  for  Eq.    (114), 
the  result  is  similar  and 

The  value  of      LA.  ^     is  determined  by  the  orthonormal  condition  im- 
posed on  the    U./yi  (Xi  *0        function.    Thus, 

|  U*  f*  2jjr  drdz     =  1. 

The  integral  is  evaluated  and  the  result  is  that 


(122) 


The  second  term  of  the  expansion  is  defined  as 
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where      cLi     =   /    KGa    UM(r,i)  Uj(rti)2jrdrdl 


and 


Xo.  =  (^A)z  +  (^/2)z  . 


The  value  of        /l^vi  ^VJ  2)     is 


CO 


/vMo;^ 


=1 


j2M 


Ac>/v>~    A  oj  J  J 


where 


cLj  -  flj  ©^©j  2tj|  */£*  tyrdrd*. 

The  solution  for       7    C  '  /  c  )         may  be  determined  by  evaluating  the 
series  of  Eq.    (121).    Thus, 


(123) 


/vi 


9^  /*» 


If  the  shape  function     JL  M    is  defined  as 


M 


TM=  ©A  + 


00 


+  •  •  I       (124; 


where 


l-*/w      '     A  oa\       /     ^^J       1  an<^       Tav     are  defined  by 


V    A/ 


Eqs.    (122),    (117),    (123),  and  (115),  the  value  of      *T  C  #,  2)       i 
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<PW     =      Bl    T^    .  W» 


The  values  of      %  ^  £)         ,  given  by  Eq .    (110),  and  Eq.    (125) 
are  used  to  evaluate  the  solvability  condition,  Eq.   (119).     Thus, 

f-^Am^M  &m    Tm    2j  dv  d*     =  O 


or 


if«t»T«rdrte    =  O. 


This  equation  is  satisfied  by  the  condition 

i  =  O  .  <125» 

When  this  result  is  substituted  into  Eq.    (118)  the  solution  for       y  (r*^) 
is  the  homogeneous  solution  to  Eq .    (118)  or  the  solution  to    Eq.    (120). 
Thus  the  value  of         ¥*(*",%)    is 


Iffr,*)      -       B/n     ±1*       . 


(127) 


The  value  of  the  perturbation  parameter,      c       /  is  chosen  to  represent 
the  average  magnitude  of  the  linear  term  in  the  expansion  of  Eq.    (113) 
weighted  by  the      JL  ^      function  and  measured  in  units  of  the  maximum 
initial  flux.    Thus 

€  ~    ;^|p  T*  [  ^%0  -  £  GO]  d*.       ws) 
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If  Eq.    (12  8)  is  differentiated  with  respect  to     £        and  the  limit  is 
taken  as     c        approaches  zero,  the  result  is 

=    /     T^     2tjt    dr  d-£ 


«M 


If  the  quantity  within  the  bracket  is  expanded  and  terms  of      &•        are 
considered  small,  then 

A«    =  /  (9^1  2rrr  dr  ds 


B 

+  a 


Av. 

oo 


+  0  (a1). 


^ 


Because  of  the  orthogonality  of  the        Tm.        function  the  second  term 
vanishes  and 


^    =      9*  ^    2jr    dr  dz 


=  1. 


Thus,      A /Y\  ~     t^m.     and  Eq.    (127)  may  be  rewritten  as 
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<P(r,i)      =A„TM. 


(129) 


30) 


The  quadratic  terms  of  the  expansion  for    TVi^A)      and     A  (£/ 
given  by  Eq.    (113)  are  evaluated  by  differentiating  Eq.    (112)  twice  with 
respect  to     t.        and  then  let     c       tend  to  zero.    The  result  is 

V*9(rtt)  +  (X  -  f  <£7?w)  <Pcw    =     (1 


<P  (r,i)     =  O  o»    iK 


The  solvability  condition  requires  that 

f   /    "      *  a  y^  ~V2  •    Z       \     -.     //  ,     (131) 

where        i    Cft^)  is  the  homogeneous  solution  to  Eq.   (130) .    This 

was  previously  determined  for  Eq.    (12  0)  and  from  that  result 

if>(r,Z)   =■   Cm    T^  ■  ("2) 


The  value  of     A       may  be  determined  by  substituting  the  correct 
functions  determined  above.    Thus, 
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X  J^  mm  X  C.m    jC^    2rjr  dr  d^      - 

The  coefficient  of      A         is  evaluated  as  follows: 

'2>  L  4sl.X£.-.  J 


V/K. 


A 


x  2tj  r  dr  d2: 


/w 


where      v-X/^        is  defined  by  Eq.    (122)  .    The  integral  on  the  right  hand 
side  is  evaluated  as  follows: 

Utfj?^  2jr  dr  d* 

7.  V 

if  terms  multiplied  by      (X         are  considered  small.    The  value  of      A 

A  =   |  yfe^  ^/©a!  2jrdr  dz.  ("si 
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For  the  fundamental  mode  (  (f\  =  X    ) 


\        - 


.561       a. 


r?      n/a7 


(134) 


The  various  terms  of  the  expansion  of  Eq.    (113)  may  be  combined 
to  give  the  equilibrium  states  of  the  reactor.    Thus,  substituting  Eqs. 
(116)  and  (12  9)  into  the  expansion  for      t  \Tt  ?j  & )       gives 


<A^o  =  A 


m 


«L+£T„  +  Od) 


(135) 


and  substituting  Eqs.    (117),    (126),  and  (133)  into  the  expansion  for 
Mw     gives 

For  the  fundamental  mode  of  the  reactor  (  //I  —  JL    )  the  integral  is 


evaluated  and 


*«>  -  C^'/rS  +tyzf    +£ 


.S63     a  - 


(136) 


+ 


It  is  apparent  that  for  any  mode  of  reactor  operation  there  will  be  two 
equilibrium  states  defined  by  Eq.    (135)  and  each  of  these  states  will  have 
unique  eigenvalues  defined  by  Eq .   (136).    However,  only  one  of  these 
values  will  represent  an  acceptable  reactor  state. 
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D.         STABILITY  ANALYSIS  OF  THE  FINITE  CYLINDRICAL  REACTOR 

The  stability  criterion  for  a  finite  cylindrical  reactor  with  negligible 

fraction  of  delayed  neutrons  is  determined  by  the  methods  previously 

considered.    An  initial  deviation,    T\       ,  due  to  temperature  feedback 

is  introduced  into  the  equilibrium  states  of  the  reactor  and  the  temporal 

effect  is  noted. 

The  governing  equation  has  the  form  given  by  Eq.    (93)  where  the 

functional  dependence  of  the  variable  is  changed  to  accommodate  the 

cylindrical  geometry 


£<fat/t) 


=  yz</W)  +[*  -  a[<P(£a  -  <£m  jfc^;o7) 


in      T> 

The  deviation  is  described  by  a  Taylor  series  expansion  about  the 
equilibrium  states,  Eq.   (135),  given  by 


(138) 


(Pfa.iSrt)  ~$(r^,e>  +7(%(w,*)  +  ^X,fr4'^^  f 


such  that 
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Eq.    (13  7)  is  differentiated  with  respect  to    ^7       and  then  the  limit  is 
taken  as  the  value  of    "?7       tends  to  zero.    The  result  of  these  operations 
yields  the  equation 

Z^Ki.t.r)     =   V^ttM,*;    +[>«)  (139) 

-  I  <ft«;\,fr,tO   -   a.  <&  c-Nj  ^("TV/O     in  J> 

A  form  for  the  solution  of     7^  (^2,f0    is  chosen  to  be 


(140) 


The  stability  of  the  system  may  be  determined  by  investigating  the  sign 
of     0\  Cw    •     If  the  sign  of    CX  (c )     were  positive,  then  the  exponential 
would  increase  with  time.    Hence,  the  contribution  of    ^Tt*       in  the 
series  of  Eq .    (138)  would  be  unbounded  with  time  and  the  system  would 
be  unstable.     On  the  other  hand,  negative  value  of    (X  (£)  would  indicate 
that  the    7^       contribution  would  die  and  stability  would  be  assured. 
The  value  of    Q\  (c/    will  be  determined  by  a  series  expansion  given  by 


The  following  will  be  directed  toward  evaluating  each  of  the  terms  of 
this  expansion. 


(141] 
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If  the  proposed  solution  for     V4   (X^t%^)    >  E<3  •   (140),  is 

substituted  intoEq.    (139)  the  result  is, 

V\x    +    [XCO    -   o(CO     +    f  <f<  174,0  d42) 

-  a    ^  <<N;J  u.  =  o   i»  s 

where    U-  ss   U.  CV  "£,  €j  •    The  first  term  of  the  series  of        &(  (£) 
is  evaluated  by  taking  the  limit  of  Eq.    (142)  as       c.      approaches  zero. 
Thus, 

V'ucy;*;     +[>0,   -   tie    -ia^T^s)]    (143) 

The  solution  to  this  equation  may  be  determined  by  the  approximation 
given  by  the  R.  Courant  and  D.  Hilbert  perturbation  technique.    A 
solution  to  Eq.    (143)  is  developed  about  the  series 

ixcr^)    =    ucr,*>    +  j  AT(^>     -h     0  (a\)  . 

The  solution  for      LA.  (  Ct  3.)    is  determined  about  the  solution  to  the 
equation 

VZU(^)      +    (^^-^o)    U(CO^    O        in  D 
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It  may  be  noted  that  the  solution  to  this  equation  is  similar  to  that  ob- 
served in  Eq.   (114).    However,  by  the  results  of  Appendix  C  one  may 
conclude  that  the  solution  for    Lk  U]  £ )  is  given  by     LK-Cf^)  =  '^yyvv  T/m 
where  the  function        //m      is  defined  in  Eq.    (115)  and       /\owc   Aom~" 0\    • 
Thus, 

Oie       -         >c^x      -      ^o/^  (144) 

and  the  linear  term  of  the      C\(€)  series  of  Eq.    (141)  is  defined.    (A 
discussion  of    (a  o      will  subsequently  be  made.)    Since  the      La(i^2) 
function  must  be  orthonormal  the  constant      <-*  /rv\     may  be  evaluated 
by  the  normality  condition 


/. 


UjL  t,     2jr    drdl    =  1  . 

The  value  of       LL  /yy\     is 


z    =  @1     = 

^      —     \J  /*V 


which  has  previously  been  determined  by  Eq .   (122).    Thus         LA  (^  "i) 
is 

The  second  term  of  the  expansion        ULOT.%)     is  defined  as 


(122) 
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oo 


f\r 


ir^)     =£ 


ci 


or 


*  Aw 


T '-T-      U.:U;a) 

A  o  aw     —    /»  o J  J 


CO 


MmCr*) 


=1 


AVt 


3- 


-V  A*. 


Aoaw    "~  No- 
where   Umj      is  defined  by  Eq.    (12  3).     The  solution  for  the       U.  (<",:£) 
function  from  Eq .   (143)  is  determined  as 

U 


^  /YH 

=      LA.  /yn 


UCr;*)    +     2  fiT^t*}  + 


oa 


] 


JL/Ws. 


(145) 


where       jL /w\     i-s  defined  by  Eq.    (124). 

The  value  of    (X  om/^     as  evaluated  from  Eq.    (144)  is  rewritten 
below 


ex  ofAftw  Ao/*\        A 


OAVL     • 


(146) 


If  the  value  of         7\0  as  given  by  Eq.    (117)  are  used,  the  result  is 


o/v\  r*\ 


For  the  fundamental  mode  (  M,  -  1    )  this  equation  reduces  to  the  for 


rn 
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uolfm  =  (if  (2J0H&  -  &)  iV$b- /m1). 

For  values  of      £       near  zero  the  values  of     OLoi/yn.    are 

0(0, »»     =0  /m   -  al 

<  -2-rjf  /m    >  A*. 

Thus  for  all  modes  of  the  deviation  except  the  fundamental  the  sign  of 

0i  0M  M\.      *s  ne<3ative  an<3  any  initial  deviation  of  the  system  will 
decay  with  time.     However,  for  the  value  /Yf\-A\    the  first  term  evalua- 
tion of  the  stability  indicates  that  there  is  a  neutral  stability  for  the 
fundamental  mode.    Thus  it  is  necessary  to  investigate  the  higher  order 
terms  of  the  expansion  of   (A  (c)    .    For  the  higher  modes  of  reactor 
operation     C\  (c  )      may  have  the  following  values, 

oi  CO      >    2tj  An  «-  m. 

Thus  the  value  of      (a(£)     may  be  positive  when  /VKL<  AV    for     /H  >  1  , 
Since  the  value  of    0(.  om  ftv\      is  so  large  when  compared  to  the  linear 
term  (order  of     £      )  of     ol  (£ )       ,  one  may  conclude  that  the  first 

term  of  the  series  is  dominant  and  the  higher  modes  are  unstable. 
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The  linear  term  in  the  expansion  of     C\  (c  J     ,  Eq.   (141),  may 
be  evaluated  by  differentiating  Eq .   (142)  with  respect  to      £       and  taking 
the  limit  as     £       tends  to  zero.    Thus, 

*  •  r^       i     ,/>  -y*      /i     (147) 

u  to*;    =   o  ck  ip  , 

The  solvability  condition  requires  that 

I  (o(  +   q&  ^C^i)  #(rtiy  ucrt2)  u(rfi)  7rjrdrd\   =0 

where      LL  Cf^)  is  the  homogeneous  solution  to  Eq .   (147).    This 

solution  was  previously  developed  for  Eq.   (143)  and  the  solution  is 

UCr,^)        =       tX/W       JL/trx    • 

.      H 

With  the  solution  for        U-  (f]^)         and  the  foiutions  for    U.ir,?)    , 

LffiCr,*)  ,  and  VVr,  tU  given  by  Eqs  .  (145),  (116),  and  (129) 
the  value  of  (A  may  be  evaluated  from  the  solvability  condition. 
Thus, 
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otU^U^J^TJL    2-rjr  dr   dz    » 

The  coefficient  of    (A        was  evaluated  for  Eq.   (128).    The  result  in- 
dicated  that  if  the  terms  of  order     (X         are  considered  very  small, 


the  integral  was  equal  to  one.     The  value  of     0\       is 

(148) 


of5?  A/j  ^m2T„  X!  2jr  dr    di. 


Since  the  fundamental  mode  of  the  reactor  is  of  greatest  interest,  the 
integral  of  Eq.    (148)  is  evaluated  for    /Wl  =  A\  sa  1    .    Thus, 


/w     -\      Q 


/  #,J^/A  2jr    dr    dH 


~     _    3,73        a_ 


The  results  of  Eqs  .    (146)  and  (148)  may  be  combined  into  the  expansion 
of        o((i)      ,  Eq.    (141), 
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<*(€.)      =      )Wm    ~    >W  (149) 


and  for  the  fundamental  mode 


From  the  previous  analysis  values  of    /Y)\   >    1.     have  been  found  to 
favor  stability.    For  the  value  of    /Yfi  -  J.    the  following  becomes 
evident, 

3.73      CL 

For  the  case  of  negative  reactivity  the  sign  of     CX  (6j     will  be  deter- 
mined by  the  value  of     C.       .A  positive  value  of      £       will  permit  a 
negative  sign  of     0\  (6)     .     The  deviation  associated  with       (X  (£ ) 
will  decrease  with  time  and  stability  results.    For  positive  reactivity 
the  sign  of  (  UL   )  is  changed  and  the  negative  value  of     <£.       permits 
stability.    The  other  states  of  reactor  equilibrium  are  unrealistic  as 
determined  by  previous  arguments  . 
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APPENDIX  A 

THE  SOLVABILITY  CONDITION  FOR  NONHOMOGENEOUS 
BOUNDARY  VALUE  PROBLEMS 

Given  the  equation 

VZ^x)     +    >)  </Vx)      -    fcx)  in  *>        (A"1} 

prove  that  this  equation  has  a  solution  if  and  only  if  the  function     TC'X) 
satisfies  the  condition: 


f  -Pm  <P<1)   dv    =  o 


where 

2    n    H  H 

V    T(x)     +     A^Ofi     ^=    O  in      D       (a-2) 


Multiplying  Eq.    (A-l)  by        i(x)        and  integrating  throughout  the 
domain  the  following  results: 

f  /        "  2  */ 

=  /    </Vx>   Ax)    dV. 
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FromEq.    (A-2)        ^7     <P(x)        ~    ""A    ^P^X)        and  this  is  substituted 
into  Eq.    (A-3) .     Thus, 

/(Wj  vVm  -  </><x;  vVrxJ   dV        <A-4> 

=  f  ¥>£)  fix)  dv  . 

By  Green's  Theorem  the  left  hand  side  of  Eq.   (A-4)  may  be  transformed 
into  a  surface  integral: 

=  /  -f(X)   <P(x)     dV. 

Thus,  if       (TiX)        is  to  be  the  solution  to  Eq.    (A-l),  the  homogeneous 
boundary  condition  on      o  J)  must  be  satisfied.    Therefore, 


/    fcx)    CPtx)     dV     «     O 


(A-5) 


Q.E.D, 
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APPENDIX  B 

A  PERTURBATION  METHOD  FOR 
LINEAR  EIGENVALUE  PROBLEMS  [7] 

The  solutions  to  equations  of  the  form 

2  -  -  (B_1) 

V  U/w  -  €  rex)  U/~    +  }\^  u^  -  O     iia  T> 

where       t     f"(X)       represents  a  small  perturbation  ,  (jc  Y'C'k)    ^    A /a  )  / 

are  provided  by  a  series  expansion  about  the  solution  of  the  equation 

V'uM      +      >/*    U/k       =    O  m    D  (B-2) 

U  /*      -    O  o*      £  D     . 

The  eigenvalues  are  held  constant  while  the  eigenfunctions  are  perturbed 
in  accordance  with  the  Taylor  series  given  by 

Um    =    U/n      +     £  AT^     +1^Wm     +•■•,  (B-3) 

If  the  series  is  substituted  into  Eq.   (B-l),  the  result  is 

-   £  rc-x)   (um  ■+■   £/\r^   +itl  laT^  +   ••) 

=   o. 
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By  separating  terms  of  the  same  power  in      £      /  the  following  equations 
are  obtained: 

Vzu„    +     AM  u^    =   O  (b-4) 

V*  at;    +    >m  atm    -    rex)  U/w       (b-5) 
V1  ur„     +     XM  urM  -    v<*)  um  .     »-6) 


Multiply  Eq.    (B-5)  by     LA  jj      and  integrate  over  the  volume,     dV     : 

/  (u*  vW  +  )\„  u^  ^)  dv  « /  ru)  u^  u*  av 


'2> 

By  Green's  Theorem  the  following  volume  integral  may  be  transformed 
into  a  surface  integral 


/     (vVm  -  at*  V'uJ  dV 


f  /I  (ujiVat^  -  /v^  VujiJ  dS  . 


However,   the  boundary  conditions  specify  that       L\o        and     /v/i\, 
are  everywhere  zero  on  the  surface  and  the  surface  integral  vanishes 
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Equation  (B-7)  may  be  rewritten  as 


I  (rf»  Vzu^    +   Ju^j/v^)  c/V    •=    d»x 

d„i    =  J  rm  u„uj,  dv. 


(B-8) 


Equation  (B-2)  reveals 

which  is  substituted  into  Eq.    (B-8)  to  give 


(B-9) 


where 


Qm^     =    /    fiTm  Ux    dV    . 


The  value  of  the  integral  Q  a>/k  is  evaluated  if  the 

normalizing  condition  /         ^m  CaV     —    X  is  imposed.     It  is 


clear  that 


/  (a/,^/v;  +  z^^-f...)2dV  «  1 .    (b-io) 

By  equating  powers  of     c.      on  the  right  and  left  sides  of  Eq.    (B-10) 

J    u^  d-x    -  1 

J   Um  /^  dy  =  O 
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or  6     &m/v\       ~     O    • 

The  function       /V/a  C  XJ        iS  determined  by  a  series  of  orthogonal 
functions  given  by 

CO 
/V^CX)       =  £      b/vnj       Uj  .  (B-ll) 

r* 

Multiplying  both  sides  by    LI  n     and  integrating  over  the  domain 

I  /V^Ujt   dV  =  /  £  b/*j  Uj  uA   dV  . 

Assuming  that  the  expansion  exists  and  that  the  integral  of  the  infinite 
series  is  equal  to  the  sum  of  the  integrals, 


e© 


a  ma    =  L  t>*^  /  uj  ujj  dv 


Ax  | 

U/V\M  =       O      . 


Since 


f  Uj  IXX  dM  =  O  j*  J! 


and 


/  uj  u^  av  -i  j=Jl    , 


D/vaa         are  evaluated  as  b  am      ==      CX,y\^    . 


the  constants         f^/         are  evaluated  as 


a        are  evaluated  as  p  /w* 


101 


I  Q  M  A 

However,  from  Eq.   (B-9)  it  is  obvious  that        D  a\  *     "~      — \ \  " 

Thus  the  series  of  orthogonal  function  given  by  Eq.  (B-ll)  may  be 
evaluated  as 


The  series  of  Eq.   (B-3)  may  be  evaluated  as 


-I.1.J-- 


J  i 


where 


clM-    =  /nxi  uw  uj  dv  . 
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APPENDIX  C 


Given  the  equation 

SJZ/i^(x)     +     >M  ^  (>r)     =   O         in    J>     (c-i) 
rf„(.x  )      =   O  ov.   «s  x> 

with  eigenfunctions  f/n  ^  ^0       anc*  eigenvalues      /\ /y\       .     Prove 

that  if 

V2cPtx)      -+•      "^^tx)       -    O  iv\D  (c-2) 


Lf>(x)        -    O  ovx     4T> 


5 


where        A    ^     A  /w      is  a  new  physical  parameter,  then  the  nontrivial 
solutions  to  Eq.    (C-2)  are 

(fix)    =   "h  ex)  >  =    > 


Av\ 


where  /Wf      and    /ft      are  independent  subscripts. 

Multiply  Eq.    (C-2)  by  IM  ^^^       and  integrate  throughout 

the  domain. 

J    fm  ( x)  V  V(x)  dV  +  >|  %  M  #*)  dV  =  O 
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By  adding  and  subtracting         ^P(X)     V     /y//»(X) 

J  (%M  V*<P(x)  -^cx)  V^<x;)   dV  (C-3) 


jf  </><*>  V2 °Lm  dV  +  XJ %b)  <P(x) 


dV    =  O  . 


From  Green's  Theorem  the  first  integral  may  be  transformed  into  a  surface 
integral  by 

f(tooV^x)-fMVztw)    dv 
■>i> 

=  [ n  ("h (■*)  V#xj  -</><*;  vf,»f>o)  ds. 


However,  from  the  homogeneous  boundary  conditions  of      %  **'  "  T  £*/  '"  ^ 
on     oJD     the  surface  integral  vanishes .    Thus  Eq.   (C-3)  is 

[<P<x>  V*^  <*>  dV  +  >  f  ^  <x>  VYx)  d  V  =  6C."4) 

FromEq.    (C-l)  V         7/VlCX)        =      -     /\/*     JmCX)     \ 

therefore 


or 


-  XM  j  <#*>  t(x^ v  +  >Y  ^  a;  <#*;  dv  -  o 

(X  -  *a0  /    ^M  ^  <*>  dV  =  O  . 


If        (A~A/rtJ^O  ,  then  the  integral  must  vanish  throughout 

the  domain  and         ^P(x)         is  orthogonal  to         /^7ii(^^  •    Hence, 
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(fix)       -     l^i  <* J 


and  the  eigenvalues  of  Eq .    (C-2)   (i.e.,      A      )  are 


Q.E.D. 
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